From 51de8202deaf132bf0d0ac56beca39fef03a4d3e Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 21 Jun 2023 06:29:36 -0400 Subject: [PATCH] +Add fluxes%tau_mag and forces%tau_mag Add new allocatable tau_mag arrays to the forcing and mech_forcing types to hold the magnitude of the wind stresses including gustiness contributions. There is also a new tau_mag diagnostic. This same information in tau_mag is being transformed into ustar, but these changes avoid division by the Boussinesq reference density (GV%Rho0), and allow for a more accurate calculation of derived fields when in non-Boussinesq mode, without having to multiply and divide by GV%Rho0. There is also a new optional tau_mag argument to extract_IOB_stresses to support these changes. These new arrays are not being used yet in the MOM6 solutions, but they are being allocated and populated in the routines that set the ustar fields, and they have been tested in changes to the modules that use ustar that will come in a subsequent commit. This commit also adds the new RLZ_T2_to_Pa element to the unit_scale_type to undo the scaling of wind stresses and it makes use of it in some of the new code. All answers are bitwise identical, but there are new arrays or elements in three transparent public types. --- .../FMS_cap/MOM_surface_forcing_gfdl.F90 | 44 ++++++++---- .../mct_cap/mom_surface_forcing_mct.F90 | 4 ++ .../nuopc_cap/mom_surface_forcing_nuopc.F90 | 4 ++ .../solo_driver/MOM_surface_forcing.F90 | 61 +++++++++++----- .../solo_driver/user_surface_forcing.F90 | 5 +- src/core/MOM.F90 | 2 + src/core/MOM_forcing_type.F90 | 70 +++++++++++++++---- src/framework/MOM_unit_scaling.F90 | 2 + 8 files changed, 144 insertions(+), 48 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 26ab6269ef..17beaecc0c 100644 --- a/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 +++ b/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 @@ -621,13 +621,16 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, endif ! Set the wind stresses and ustar. - if (associated(fluxes%ustar) .and. associated(fluxes%ustar_gustless)) then + if (associated(fluxes%ustar) .and. associated(fluxes%ustar_gustless) .and. associated(fluxes%tau_mag)) then call extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, ustar=fluxes%ustar, & - gustless_ustar=fluxes%ustar_gustless) - elseif (associated(fluxes%ustar)) then - call extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, ustar=fluxes%ustar) - elseif (associated(fluxes%ustar_gustless)) then - call extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, gustless_ustar=fluxes%ustar_gustless) + mag_tau=fluxes%tau_mag, gustless_ustar=fluxes%ustar_gustless) + else + if (associated(fluxes%ustar)) & + call extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, ustar=fluxes%ustar) + if (associated(fluxes%ustar_gustless)) & + call extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, gustless_ustar=fluxes%ustar_gustless) + if (associated(fluxes%tau_mag)) & + call extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, mag_tau=fluxes%tau_mag) endif if (coupler_type_initialized(fluxes%tr_fluxes) .and. & @@ -671,7 +674,8 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_ real, dimension(SZI_(G),SZJ_(G)) :: & 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]. + 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] 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] @@ -775,12 +779,13 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_ ! Set the wind stresses and ustar. if (wt1 <= 0.0) then call extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux=forces%taux, tauy=forces%tauy, & - ustar=forces%ustar, tau_halo=1) + ustar=forces%ustar, mag_tau=forces%tau_mag, tau_halo=1) else call extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux=forces%taux, tauy=forces%tauy, & - ustar=ustar_tmp, tau_halo=1) + ustar=ustar_tmp, mag_tau=tau_mag_tmp, tau_halo=1) do j=js,je ; do i=is,ie forces%ustar(i,j) = wt1*forces%ustar(i,j) + wt2*ustar_tmp(i,j) + forces%tau_mag(i,j) = wt1*forces%tau_mag(i,j) + wt2*tau_mag_tmp(i,j) enddo ; enddo endif @@ -877,7 +882,7 @@ end subroutine convert_IOB_to_forces !! Ice_ocean_boundary_type into optional argument arrays, including changes of units, sign !! conventions, and putting the fields into arrays with MOM-standard sized halos. subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy, ustar, & - gustless_ustar, tau_halo) + gustless_ustar, mag_tau, tau_halo) type(ice_ocean_boundary_type), & target, intent(in) :: IOB !< An ice-ocean boundary type with fluxes to drive !! the ocean in a coupled model @@ -897,6 +902,9 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy, real, dimension(SZI_(G),SZJ_(G)), & optional, intent(out) :: gustless_ustar !< The surface friction velocity without !! 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] integer, optional, intent(in) :: tau_halo !< The halo size of wind stresses to set, 0 by default. ! Local variables @@ -913,7 +921,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy, real :: tau_mag ! magnitude of the wind stress [R Z L T-2 ~> Pa] real :: stress_conversion ! A unit conversion factor from Pa times any stress multiplier [R Z L T-2 Pa-1 ~> 1] - logical :: do_ustar, do_gustless + logical :: do_ustar, do_gustless, do_tau_mag integer :: wind_stagger ! AGRID, BGRID_NE, or CGRID_NE (integers from MOM_domains) integer :: i, j, is, ie, js, je, ish, ieh, jsh, jeh, Isqh, Ieqh, Jsqh, Jeqh, i0, j0, halo @@ -926,7 +934,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy, IRho0 = US%L_to_Z / CS%Rho0 stress_conversion = US%Pa_to_RLZ_T2 * CS%wind_stress_multiplier - do_ustar = present(ustar) ; do_gustless = present(gustless_ustar) + do_ustar = present(ustar) ; do_gustless = present(gustless_ustar) ; do_tau_mag = present(mag_tau) wind_stagger = CS%wind_stagger if ((IOB%wind_stagger == AGRID) .or. (IOB%wind_stagger == BGRID_NE) .or. & @@ -1019,13 +1027,13 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy, endif ! endif for extracting wind stress fields with various staggerings endif - if (do_ustar .or. do_gustless) then + if (do_ustar .or. do_tau_mag .or. do_gustless) then ! Set surface friction velocity directly or as a function of staggering. ! ustar is required for the bulk mixed layer formulation and other turbulent mixing ! 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 - if (do_ustar) then ; do j=js,je ; do i=is,ie + 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 if ((wind_stagger == CGRID_NE) .or. & @@ -1035,7 +1043,10 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy, (G%mask2dBu(I,J-1) + G%mask2dBu(I-1,J))) > 0.0)) ) & gustiness = CS%gust(i,j) endif - ustar(i,j) = sqrt(gustiness*IRho0 + IRho0*US%Pa_to_RLZ_T2*IOB%stress_mag(i-i0,j-j0)) + if (do_tau_mag) & + mag_tau(i,j) = gustiness + US%Pa_to_RLZ_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)) enddo ; enddo ; endif if (CS%answer_date < 20190101) then if (do_gustless) then ; do j=js,je ; do i=is,ie @@ -1059,6 +1070,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy, if (CS%read_gust_2d) gustiness = CS%gust(i,j) endif if (do_ustar) ustar(i,j) = sqrt(gustiness*IRho0 + IRho0 * tau_mag) + if (do_tau_mag) mag_tau(i,j) = gustiness + tau_mag if (CS%answer_date < 20190101) then if (do_gustless) gustless_ustar(i,j) = sqrt(US%L_to_Z*tau_mag / CS%Rho0) else @@ -1071,6 +1083,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy, 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 (CS%answer_date < 20190101) then if (do_gustless) gustless_ustar(i,j) = sqrt(US%L_to_Z*tau_mag / CS%Rho0) else @@ -1092,6 +1105,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy, if (CS%read_gust_2d) 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 (CS%answer_date < 20190101) then if (do_gustless) gustless_ustar(i,j) = sqrt(US%L_to_Z*tau_mag / CS%Rho0) else diff --git a/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 b/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 index 0364d46ddc..0ba3d5242b 100644 --- a/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 +++ b/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 @@ -771,6 +771,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%ustar(i,j) = sqrt(gustiness*Irho0 + Irho0*tau_mag) enddo ; enddo @@ -796,6 +797,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%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 @@ -817,8 +819,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%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%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 2c8e3db8bd..f774ae91a8 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -845,6 +845,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%ustar(i,j) = sqrt(gustiness*Irho0 + Irho0*tau_mag) enddo ; enddo call pass_vector(forces%taux, forces%tauy, G%Domain, halo=1) @@ -870,6 +871,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%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 @@ -891,8 +893,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%ustar(i,j) = sqrt(CS%gust(i,j)*Irho0 + Irho0*sqrt(taux2 + tauy2)) else + forces%tau_mag(i,j) = CS%gust(i,j) + 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 0e8aedb8d0..859bfd81c8 100644 --- a/config_src/drivers/solo_driver/MOM_surface_forcing.F90 +++ b/config_src/drivers/solo_driver/MOM_surface_forcing.F90 @@ -407,10 +407,16 @@ subroutine wind_forcing_const(sfc_state, forces, tau_x0, tau_y0, day, G, US, CS) 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 ) 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 ) 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 + enddo ; enddo ; endif endif call callTree_leave("wind_forcing_const") @@ -522,9 +528,11 @@ subroutine wind_forcing_gyres(sfc_state, forces, day, G, US, CS) ! set the friction velocity if (CS%answer_date < 20190101) 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%ustar(i,j) = sqrt(US%L_to_Z * ((CS%gust_const/CS%Rho0) + & - sqrt(0.5*(forces%tauy(i,j-1)*forces%tauy(i,j-1) + forces%tauy(i,j)*forces%tauy(i,j) + & - forces%taux(i-1,j)*forces%taux(i-1,j) + forces%taux(i,j)*forces%taux(i,j)))/CS%Rho0) ) + sqrt(0.5*(forces%tauy(i,J-1)*forces%tauy(i,J-1) + forces%tauy(i,J)*forces%tauy(i,J) + & + forces%taux(I-1,j)*forces%taux(I-1,j) + forces%taux(I,j)*forces%taux(I,j)))/CS%Rho0) ) enddo ; enddo else call stresses_to_ustar(forces, G, US, CS) @@ -725,11 +733,12 @@ subroutine wind_forcing_from_file(sfc_state, forces, day, G, US, CS) if (.not.read_Ustar) then if (CS%read_gust_2d) then do j=js,je ; do i=is,ie - forces%ustar(i,j) = sqrt((CS%gust(i,j) + & - sqrt(temp_x(i,j)*temp_x(i,j) + temp_y(i,j)*temp_y(i,j))) * US%L_to_Z / CS%Rho0) + forces%tau_mag(i,j) = CS%gust(i,j) + sqrt(temp_x(i,j)**2 + temp_y(i,j)**2) + forces%ustar(i,j) = sqrt(forces%tau_mag(i,j) * US%L_to_Z / CS%Rho0) enddo ; enddo else 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%ustar(i,j) = sqrt(US%L_to_Z * (CS%gust_const/CS%Rho0 + & sqrt(temp_x(i,j)*temp_x(i,j) + temp_y(i,j)*temp_y(i,j)) / CS%Rho0) ) enddo ; enddo @@ -772,15 +781,19 @@ subroutine wind_forcing_from_file(sfc_state, forces, day, G, US, CS) if (.not.read_Ustar) then if (CS%read_gust_2d) then do j=js,je ; do i=is,ie - forces%ustar(i,j) = sqrt((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 / CS%Rho0 ) + 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))) + forces%ustar(i,j) = sqrt( forces%tau_mag(i,j) * US%L_to_Z / CS%Rho0 ) enddo ; enddo else 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%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)) + 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 @@ -792,6 +805,9 @@ subroutine wind_forcing_from_file(sfc_state, forces, day, G, US, CS) if (read_Ustar) then call MOM_read_data(filename, CS%Ustar_var, forces%ustar(:,:), & G%Domain, timelevel=time_lev, scale=US%m_to_Z*US%T_to_s) + do j=js,je ; do i=is,ie + forces%tau_mag(i,j) = US%Z_to_L * CS%Rho0 * forces%ustar(i,j)**2 + enddo ; enddo endif CS%wind_last_lev = time_lev @@ -815,6 +831,7 @@ subroutine wind_forcing_by_data_override(sfc_state, forces, day, G, US, CS) ! Local variables real :: temp_x(SZI_(G),SZJ_(G)) ! Pseudo-zonal wind stresses at h-points [R Z L T-2 ~> Pa]. real :: temp_y(SZI_(G),SZJ_(G)) ! Pseudo-meridional wind stresses at h-points [R Z L T-2 ~> Pa]. + real :: ustar_tmp(SZI_(G),SZJ_(G)) ! The pre-override value of ustar [Z T-1 ~> m s-1] integer :: i, j call callTree_enter("wind_forcing_by_data_override, MOM_surface_forcing.F90") @@ -840,17 +857,25 @@ subroutine wind_forcing_by_data_override(sfc_state, forces, day, G, US, CS) if (CS%read_gust_2d) then call data_override(G%Domain, 'gust', CS%gust, day, scale=US%Pa_to_RLZ_T2) do j=G%jsc,G%jec ; do i=G%isc,G%iec - forces%ustar(i,j) = sqrt((sqrt(temp_x(i,j)**2 + temp_y(i,j)**2) + & - CS%gust(i,j)) * US%L_to_Z / CS%Rho0) + forces%tau_mag(i,j) = sqrt(temp_x(i,j)**2 + temp_y(i,j)**2) + CS%gust(i,j) + forces%ustar(i,j) = sqrt( forces%tau_mag(i,j) * US%L_to_Z / CS%Rho0 ) enddo ; enddo else 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 + ! forces%ustar(i,j) = sqrt( forces%tau_mag(i,j) * US%L_to_Z / CS%Rho0 ) forces%ustar(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 + ! Give the data override the option to modify the newly calculated forces%ustar. + ustar_tmp(:,:) = forces%ustar(:,:) call data_override(G%Domain, 'ustar', forces%ustar, day, scale=US%m_to_Z*US%T_to_s) + ! Only reset values where data override of ustar has occurred + do j=G%jsc,G%jec ; do i=G%isc,G%iec ; if (ustar_tmp(i,j) /= forces%ustar(i,j)) then + forces%tau_mag(i,j) = US%Z_to_L * CS%Rho0 * forces%ustar(i,j)**2 + endif ; enddo ; enddo call pass_vector(forces%taux, forces%tauy, G%Domain, To_All) @@ -875,15 +900,17 @@ subroutine stresses_to_ustar(forces, G, US, CS) if (CS%read_gust_2d) then do j=js,je ; do i=is,ie - forces%ustar(i,j) = sqrt( (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))) ) * I_rho ) + 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))) + forces%ustar(i,j) = sqrt( forces%tau_mag(i,j) * I_rho ) enddo ; enddo else do j=js,je ; do i=is,ie - forces%ustar(i,j) = sqrt( (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))) ) * I_rho ) + 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%ustar(i,j) = sqrt( forces%tau_mag(i,j) * I_rho ) enddo ; enddo endif diff --git a/config_src/drivers/solo_driver/user_surface_forcing.F90 b/config_src/drivers/solo_driver/user_surface_forcing.F90 index 42e732bb73..d7d3b89a8a 100644 --- a/config_src/drivers/solo_driver/user_surface_forcing.F90 +++ b/config_src/drivers/solo_driver/user_surface_forcing.F90 @@ -88,9 +88,10 @@ subroutine USER_wind_forcing(sfc_state, forces, day, G, US, CS) ! is always positive. 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((CS%gust_const + & + 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/CS%Rho0)) + 0.5*(forces%tauy(i,J-1)**2 + forces%tauy(i,J)**2))) + forces%ustar(i,j) = G%mask2dT(i,j) * sqrt(forces%tau_mag(i,j) * (US%L_to_Z/CS%Rho0)) enddo ; enddo ; endif end subroutine USER_wind_forcing diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index af8481fd1c..7b9f2f9d3f 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -651,6 +651,8 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS call create_group_pass(pass_tau_ustar_psurf, forces%taux, forces%tauy, G%Domain) if (associated(forces%ustar)) & call create_group_pass(pass_tau_ustar_psurf, forces%ustar, G%Domain) + if (associated(forces%tau_mag)) & + call create_group_pass(pass_tau_ustar_psurf, forces%tau_mag, G%Domain) if (associated(forces%p_surf)) & call create_group_pass(pass_tau_ustar_psurf, forces%p_surf, G%Domain) if (G%nonblocking_updates) then diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 9a5e1f48f5..a59c33d525 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -68,6 +68,9 @@ module MOM_forcing_type ! surface stress components and turbulent velocity scale real, pointer, dimension(:,:) :: & 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] ustar_gustless => NULL() !< surface friction velocity scale without any !! any augmentation for gustiness [Z T-1 ~> m s-1]. @@ -222,6 +225,8 @@ module MOM_forcing_type real, pointer, dimension(:,:) :: & taux => NULL(), & !< zonal wind stress [R L Z T-2 ~> Pa] tauy => NULL(), & !< meridional wind stress [R L Z T-2 ~> Pa] + 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] ustar => NULL(), & !< surface friction velocity scale [Z T-1 ~> m s-1]. net_mass_src => NULL() !< The net mass source to the ocean [R Z T-1 ~> kg m-2 s-1] @@ -359,6 +364,7 @@ module MOM_forcing_type integer :: id_taux = -1 integer :: id_tauy = -1 integer :: id_ustar = -1 + integer :: id_tau_mag = -1 integer :: id_psurf = -1 integer :: id_TKE_tidal = -1 @@ -1083,6 +1089,8 @@ subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift) ! and js...je as their extent. if (associated(fluxes%ustar)) & call hchksum(fluxes%ustar, mesg//" fluxes%ustar", G%HI, haloshift=hshift, scale=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, scale=US%RLZ_T2_to_Pa) if (associated(fluxes%buoy)) & call hchksum(fluxes%buoy, mesg//" fluxes%buoy ", G%HI, haloshift=hshift, scale=US%L_to_m**2*US%s_to_T**3) if (associated(fluxes%sw)) & @@ -1186,11 +1194,13 @@ subroutine MOM_mech_forcing_chksum(mesg, forces, G, US, haloshift) ! and js...je as their extent. if (associated(forces%taux) .and. associated(forces%tauy)) & call uvchksum(mesg//" forces%tau[xy]", forces%taux, forces%tauy, G%HI, & - haloshift=hshift, symmetric=.true., scale=US%RZ_T_to_kg_m2s*US%L_T_to_m_s) + haloshift=hshift, symmetric=.true., scale=US%RLZ_T2_to_Pa) if (associated(forces%p_surf)) & call hchksum(forces%p_surf, mesg//" forces%p_surf", G%HI, haloshift=hshift, scale=US%RL2_T2_to_Pa) if (associated(forces%ustar)) & call hchksum(forces%ustar, mesg//" forces%ustar", G%HI, haloshift=hshift, scale=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, scale=US%RLZ_T2_to_Pa) 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., & @@ -1237,6 +1247,7 @@ subroutine forcing_SinglePointPrint(fluxes, G, i, j, mesg) write(0,'(2a)') 'MOM_forcing_type, forcing_SinglePointPrint: Called from ',mesg write(0,'(a,2es15.3)') 'MOM_forcing_type, forcing_SinglePointPrint: lon,lat = ',G%geoLonT(i,j),G%geoLatT(i,j) call locMsg(fluxes%ustar,'ustar') + call locMsg(fluxes%tau_mag,'tau_mag') call locMsg(fluxes%buoy,'buoy') call locMsg(fluxes%sw,'sw') call locMsg(fluxes%sw_vis_dir,'sw_vis_dir') @@ -1305,18 +1316,22 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, handles%id_taux = register_diag_field('ocean_model', 'taux', diag%axesCu1, Time, & 'Zonal surface stress from ocean interactions with atmos and ice', & - 'Pa', conversion=US%RZ_T_to_kg_m2s*US%L_T_to_m_s, & + 'Pa', conversion=US%RLZ_T2_to_Pa, & standard_name='surface_downward_x_stress', cmor_field_name='tauuo', & cmor_units='N m-2', cmor_long_name='Surface Downward X Stress', & cmor_standard_name='surface_downward_x_stress') handles%id_tauy = register_diag_field('ocean_model', 'tauy', diag%axesCv1, Time, & 'Meridional surface stress ocean interactions with atmos and ice', & - 'Pa', conversion=US%RZ_T_to_kg_m2s*US%L_T_to_m_s, & + 'Pa', conversion=US%RLZ_T2_to_Pa, & standard_name='surface_downward_y_stress', cmor_field_name='tauvo', & cmor_units='N m-2', cmor_long_name='Surface Downward Y Stress', & cmor_standard_name='surface_downward_y_stress') + 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) + handles%id_ustar = register_diag_field('ocean_model', 'ustar', diag%axesT1, Time, & 'Surface friction velocity = [(gustiness + tau_magnitude)/rho0]^(1/2)', & 'm s-1', conversion=US%Z_to_m*US%s_to_T) @@ -2046,6 +2061,7 @@ subroutine fluxes_accumulate(flux_tmp, fluxes, G, wt2, forces) fluxes%p_surf_full(i,j) = forces%p_surf_full(i,j) fluxes%ustar(i,j) = wt1*fluxes%ustar(i,j) + wt2*forces%ustar(i,j) + fluxes%tau_mag(i,j) = wt1*fluxes%tau_mag(i,j) + wt2*forces%tau_mag(i,j) enddo ; enddo else do j=js,je ; do i=is,ie @@ -2053,6 +2069,7 @@ subroutine fluxes_accumulate(flux_tmp, fluxes, G, wt2, forces) fluxes%p_surf_full(i,j) = flux_tmp%p_surf_full(i,j) fluxes%ustar(i,j) = wt1*fluxes%ustar(i,j) + wt2*flux_tmp%ustar(i,j) + fluxes%tau_mag(i,j) = wt1*fluxes%tau_mag(i,j) + wt2*flux_tmp%tau_mag(i,j) enddo ; enddo endif @@ -2173,6 +2190,12 @@ subroutine copy_common_forcing_fields(forces, fluxes, G, skip_pres) enddo ; enddo endif + if (associated(forces%tau_mag) .and. associated(fluxes%tau_mag)) then + do j=js,je ; do i=is,ie + fluxes%tau_mag(i,j) = forces%tau_mag(i,j) + enddo ; enddo + endif + if (do_pres) then if (associated(forces%p_surf) .and. associated(fluxes%p_surf)) then do j=js,je ; do i=is,ie @@ -2304,6 +2327,12 @@ subroutine copy_back_forcing_fields(fluxes, forces, G) enddo ; enddo endif + if (associated(forces%tau_mag) .and. associated(fluxes%tau_mag)) then + do j=js,je ; do i=is,ie + forces%tau_mag(i,j) = fluxes%tau_mag(i,j) + enddo ; enddo + endif + end subroutine copy_back_forcing_fields !> Offer mechanical forcing fields for diagnostics for those @@ -2946,6 +2975,9 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h if ((handles%id_buoy > 0) .and. associated(fluxes%buoy)) & call post_data(handles%id_buoy, fluxes%buoy, diag) + if ((handles%id_tau_mag > 0) .and. associated(fluxes%tau_mag)) & + call post_data(handles%id_tau_mag, fluxes%tau_mag, diag) + if ((handles%id_ustar > 0) .and. associated(fluxes%ustar)) & call post_data(handles%id_ustar, fluxes%ustar, diag) @@ -3015,6 +3047,7 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & call myAlloc(fluxes%ustar,isd,ied,jsd,jed, ustar) call myAlloc(fluxes%ustar_gustless,isd,ied,jsd,jed, ustar) + call myAlloc(fluxes%tau_mag,isd,ied,jsd,jed, ustar) call myAlloc(fluxes%evap,isd,ied,jsd,jed, water) call myAlloc(fluxes%lprec,isd,ied,jsd,jed, water) @@ -3150,6 +3183,7 @@ subroutine allocate_mech_forcing_by_group(G, forces, stress, ustar, shelf, & call myAlloc(forces%tauy,isd,ied,JsdB,JedB, stress) call myAlloc(forces%ustar,isd,ied,jsd,jed, ustar) + call myAlloc(forces%tau_mag,isd,ied,jsd,jed, ustar) call myAlloc(forces%p_surf,isd,ied,jsd,jed, press) call myAlloc(forces%p_surf_full,isd,ied,jsd,jed, press) @@ -3218,8 +3252,7 @@ subroutine get_forcing_groups(fluxes, water, heat, ustar, press, shelf, & ! to some degree. But since this would be enforced at the driver level, ! we handle them here as independent flags. - ustar = associated(fluxes%ustar) & - .and. associated(fluxes%ustar_gustless) + ustar = associated(fluxes%ustar) .and. associated(fluxes%ustar_gustless) ! TODO: Check for all associated fields, but for now just check one as a marker water = associated(fluxes%evap) heat = associated(fluxes%seaice_melt_heat) @@ -3276,6 +3309,7 @@ subroutine deallocate_forcing_type(fluxes) if (associated(fluxes%ustar)) deallocate(fluxes%ustar) if (associated(fluxes%ustar_gustless)) deallocate(fluxes%ustar_gustless) + if (associated(fluxes%tau_mag)) deallocate(fluxes%tau_mag) if (associated(fluxes%buoy)) deallocate(fluxes%buoy) if (associated(fluxes%sw)) deallocate(fluxes%sw) if (associated(fluxes%seaice_melt_heat)) deallocate(fluxes%seaice_melt_heat) @@ -3334,9 +3368,10 @@ end subroutine deallocate_forcing_type subroutine deallocate_mech_forcing(forces) type(mech_forcing), intent(inout) :: forces !< Forcing fields structure - if (associated(forces%taux)) deallocate(forces%taux) - if (associated(forces%tauy)) deallocate(forces%tauy) - if (associated(forces%ustar)) deallocate(forces%ustar) + if (associated(forces%taux)) deallocate(forces%taux) + if (associated(forces%tauy)) deallocate(forces%tauy) + if (associated(forces%ustar)) deallocate(forces%ustar) + if (associated(forces%tau_mag)) deallocate(forces%tau_mag) if (associated(forces%p_surf)) deallocate(forces%p_surf) if (associated(forces%p_surf_full)) deallocate(forces%p_surf_full) if (associated(forces%net_mass_src)) deallocate(forces%net_mass_src) @@ -3365,6 +3400,7 @@ subroutine rotate_forcing(fluxes_in, fluxes, turns) if (do_ustar) then call rotate_array(fluxes_in%ustar, turns, fluxes%ustar) call rotate_array(fluxes_in%ustar_gustless, turns, fluxes%ustar_gustless) + call rotate_array(fluxes_in%tau_mag, turns, fluxes%tau_mag) endif if (do_water) then @@ -3495,8 +3531,10 @@ subroutine rotate_mech_forcing(forces_in, turns, forces) call rotate_vector(forces_in%taux, forces_in%tauy, turns, & forces%taux, forces%tauy) - if (do_ustar) & + if (do_ustar) then call rotate_array(forces_in%ustar, turns, forces%ustar) + call rotate_array(forces_in%tau_mag, turns, forces%tau_mag) + endif if (do_shelf) then call rotate_array_pair( & @@ -3555,24 +3593,27 @@ subroutine homogenize_mech_forcing(forces, G, US, Rho0, UpdateUstar) do_press, do_iceberg) if (do_stress) then - tx_mean = global_area_mean_u(forces%taux, G, tmp_scale=US%Z_to_L*US%RL2_T2_to_Pa) + tx_mean = global_area_mean_u(forces%taux, G, tmp_scale=US%RLZ_T2_to_Pa) do j=js,je ; do i=isB,ieB if (G%mask2dCu(I,j) > 0.0) forces%taux(I,j) = tx_mean enddo ; enddo - ty_mean = global_area_mean_v(forces%tauy, G, tmp_scale=US%Z_to_L*US%RL2_T2_to_Pa) + ty_mean = global_area_mean_v(forces%tauy, G, tmp_scale=US%RLZ_T2_to_Pa) do j=jsB,jeB ; do i=is,ie if (G%mask2dCv(i,J) > 0.0) forces%tauy(i,J) = ty_mean enddo ; enddo if (tau2ustar) then - do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) > 0.0) forces%ustar(i,j) = sqrt(sqrt(tx_mean**2 + ty_mean**2)*Irho0) - enddo ; enddo + do j=js,je ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then + forces%tau_mag(i,j) = sqrt(tx_mean**2 + ty_mean**2) + forces%ustar(i,j) = sqrt(forces%tau_mag(i,j) * Irho0) + endif ; enddo ; enddo else call homogenize_field_t(forces%ustar, G, tmp_scale=US%Z_to_m*US%s_to_T) + call homogenize_field_t(forces%tau_mag, G, tmp_scale=US%RLZ_T2_to_Pa) endif else if (do_ustar) then call homogenize_field_t(forces%ustar, G, tmp_scale=US%Z_to_m*US%s_to_T) + call homogenize_field_t(forces%tau_mag, G, tmp_scale=US%RLZ_T2_to_Pa) endif endif @@ -3613,6 +3654,7 @@ subroutine homogenize_forcing(fluxes, G, GV, US) if (do_ustar) then call homogenize_field_t(fluxes%ustar, G, tmp_scale=US%Z_to_m*US%s_to_T) call homogenize_field_t(fluxes%ustar_gustless, G, tmp_scale=US%Z_to_m*US%s_to_T) + call homogenize_field_t(fluxes%tau_mag, G, tmp_scale=US%RLZ_T2_to_Pa) endif if (do_water) then diff --git a/src/framework/MOM_unit_scaling.F90 b/src/framework/MOM_unit_scaling.F90 index 482c2eec7a..868352102e 100644 --- a/src/framework/MOM_unit_scaling.F90 +++ b/src/framework/MOM_unit_scaling.F90 @@ -52,6 +52,7 @@ module MOM_unit_scaling 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] real :: W_m2_to_RZ3_T3 !< Convert turbulent kinetic energy fluxes from W m-2 to R Z3 T-3 [R Z3 m2 T-3 W-1 ~> 1] real :: RL2_T2_to_Pa !< Convert pressures from R L2 T-2 to Pa [Pa T2 R-1 L-2 ~> 1] + real :: RLZ_T2_to_Pa !< Convert wind stresses from R L Z T-2 to Pa [Pa T2 R-1 L-1 Z-1 ~> 1] real :: Pa_to_RL2_T2 !< Convert pressures from Pa to R L2 T-2 [R L2 T-2 Pa-1 ~> 1] real :: Pa_to_RLZ_T2 !< Convert wind stresses from Pa to R L Z T-2 [R L Z T-2 Pa-1 ~> 1] @@ -221,6 +222,7 @@ subroutine set_unit_scaling_combos(US) US%RL2_T2_to_Pa = US%R_to_kg_m3 * US%L_T_to_m_s**2 US%Pa_to_RL2_T2 = US%kg_m3_to_R * US%m_s_to_L_T**2 ! 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 end subroutine set_unit_scaling_combos