Skip to content

Commit

Permalink
+Rescaled units of forces%taux to [R Z L T-2]
Browse files Browse the repository at this point in the history 
  Rescaled the units for the wind stresses, forces%taux and fluxes%tauy, from
[Pa] to [R Z L T-2], for expanded dimensional consistency testing.  All answers
are bitwise identical, but there are changes in the dimensions of two elements
in a transparent public type.  Some changes in the mct_driver and the
nuopc_driver are not well tested, but are analogous to changes in well-tested
code.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Oct 3, 2019
1 parent a669294 commit 1cbf498
Show file tree
Hide file tree
Showing 14 changed files with 179 additions and 135 deletions.
61 changes: 33 additions & 28 deletions config_src/coupled_driver/MOM_surface_forcing_gfdl.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -862,9 +862,9 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
type(surface_forcing_CS),pointer :: CS !< A pointer to the control structure returned by a
!! previous call to surface_forcing_init.
real, dimension(SZIB_(G),SZJ_(G)), &
optional, intent(inout) :: taux !< The zonal wind stresses on a C-grid [Pa].
optional, intent(inout) :: taux !< The zonal wind stresses on a C-grid [R Z L T-2 ~> Pa].
real, dimension(SZI_(G),SZJB_(G)), &
optional, intent(inout) :: tauy !< The meridional wind stresses on a C-grid [Pa].
optional, intent(inout) :: tauy !< The meridional wind stresses on a C-grid [R Z L T-2 ~> Pa].
real, dimension(SZI_(G),SZJ_(G)), &
optional, intent(inout) :: ustar !< The surface friction velocity [Z T-1 ~> m s-1].
real, dimension(SZI_(G),SZJ_(G)), &
Expand All @@ -873,17 +873,18 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
integer, optional, intent(in) :: tau_halo !< The halo size of wind stresses to set, 0 by default.

! Local variables
real, dimension(SZI_(G),SZJ_(G)) :: taux_in_A ! Zonal wind stresses [Pa] at h points
real, dimension(SZI_(G),SZJ_(G)) :: tauy_in_A ! Meridional wind stresses [Pa] at h points
real, dimension(SZI_(G),SZJ_(G)) :: taux_in_A ! Zonal wind stresses [R Z L T-2 ~> Pa] at h points
real, dimension(SZI_(G),SZJ_(G)) :: tauy_in_A ! Meridional wind stresses [R Z L T-2 ~> Pa] at h points
real, dimension(SZIB_(G),SZJ_(G)) :: taux_in_C ! Zonal wind stresses [Pa] at u points
real, dimension(SZI_(G),SZJB_(G)) :: tauy_in_C ! Meridional wind stresses [Pa] at v points
real, dimension(SZI_(G),SZJB_(G)) :: tauy_in_C ! Meridional wind stresses [R Z L T-2 ~> Pa] at v points
real, dimension(SZIB_(G),SZJB_(G)) :: taux_in_B ! Zonal wind stresses [Pa] at q points
real, dimension(SZIB_(G),SZJB_(G)) :: tauy_in_B ! Meridional wind stresses [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 [Pa]
real :: gustiness ! unresolved gustiness that contributes to ustar [R Z L T-2 ~> Pa]
real :: Irho0 ! Inverse of the mean density rescaled to [Z2 s2 m T-2 kg-1 ~> m3 kg-1]
real :: taux2, tauy2 ! squared wind stresses [Pa2]
real :: tau_mag ! magnitude of the wind stress [Pa]
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 :: 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
integer :: wind_stagger ! AGRID, BGRID_NE, or CGRID_NE (integers from MOM_domains)
Expand All @@ -896,6 +897,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
i0 = is - index_bounds(1) ; j0 = js - index_bounds(3)

Irho0 = (US%m_to_Z*US%T_to_s)**2 / CS%Rho0
stress_conversion = US%kg_m3_to_R*US%m_s_to_L_T**2*US%L_to_Z * CS%wind_stress_multiplier

do_ustar = present(ustar) ; do_gustless = present(gustless_ustar)

Expand All @@ -916,8 +918,8 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
taux_in_B(:,:) = 0.0 ; tauy_in_B(:,:) = 0.0
if (associated(IOB%u_flux).and.associated(IOB%v_flux)) then
do J=js,je ; do I=is,ie
taux_in_B(I,J) = IOB%u_flux(i-i0,j-j0) * CS%wind_stress_multiplier
tauy_in_B(I,J) = IOB%v_flux(i-i0,j-j0) * CS%wind_stress_multiplier
taux_in_B(I,J) = IOB%u_flux(i-i0,j-j0) * stress_conversion
tauy_in_B(I,J) = IOB%v_flux(i-i0,j-j0) * stress_conversion
enddo ; enddo
endif

Expand All @@ -942,8 +944,8 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
taux_in_A(:,:) = 0.0 ; tauy_in_A(:,:) = 0.0
if (associated(IOB%u_flux).and.associated(IOB%v_flux)) then
do j=js,je ; do i=is,ie
taux_in_A(i,j) = IOB%u_flux(i-i0,j-j0) * CS%wind_stress_multiplier
tauy_in_A(i,j) = IOB%v_flux(i-i0,j-j0) * CS%wind_stress_multiplier
taux_in_A(i,j) = IOB%u_flux(i-i0,j-j0) * stress_conversion
tauy_in_A(i,j) = IOB%v_flux(i-i0,j-j0) * stress_conversion
enddo ; enddo
endif

Expand Down Expand Up @@ -971,8 +973,8 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
taux_in_C(:,:) = 0.0 ; tauy_in_C(:,:) = 0.0
if (associated(IOB%u_flux).and.associated(IOB%v_flux)) then
do j=js,je ; do i=is,ie
taux_in_C(I,j) = IOB%u_flux(i-i0,j-j0) * CS%wind_stress_multiplier
tauy_in_C(i,J) = IOB%v_flux(i-i0,j-j0) * CS%wind_stress_multiplier
taux_in_C(I,j) = IOB%u_flux(i-i0,j-j0) * stress_conversion
tauy_in_C(i,J) = IOB%v_flux(i-i0,j-j0) * stress_conversion
enddo ; enddo
endif

Expand Down Expand Up @@ -1029,23 +1031,23 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
((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
if (do_ustar) ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * tau_mag)
if (do_ustar) ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * US%R_to_kg_m3*US%L_T_to_m_s**2*US%Z_to_L*tau_mag)
if (CS%answers_2018) then
if (do_gustless) gustless_ustar(i,j) = US%m_to_Z*US%T_to_s * sqrt(tau_mag / CS%Rho0)
if (do_gustless) gustless_ustar(i,j) = sqrt(US%R_to_kg_m3*US%L_to_Z*tau_mag / CS%Rho0)
else
if (do_gustless) gustless_ustar(i,j) = sqrt(Irho0 * tau_mag)
if (do_gustless) gustless_ustar(i,j) = sqrt(Irho0 * US%R_to_kg_m3*US%L_T_to_m_s**2*US%Z_to_L*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)
gustiness = CS%gust_const
if (CS%read_gust_2d .and. (G%mask2dT(i,j) > 0)) gustiness = CS%gust(i,j)
if (do_ustar) ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * tau_mag)
if (do_ustar) ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * US%R_to_kg_m3*US%L_T_to_m_s**2*US%Z_to_L*tau_mag)
if (CS%answers_2018) then
if (do_gustless) gustless_ustar(i,j) = US%m_to_Z*US%T_to_s * sqrt(tau_mag / CS%Rho0)
if (do_gustless) gustless_ustar(i,j) = sqrt(US%R_to_kg_m3*US%L_to_Z*tau_mag / CS%Rho0)
else
if (do_gustless) gustless_ustar(i,j) = sqrt(Irho0 * tau_mag)
if (do_gustless) gustless_ustar(i,j) = sqrt(Irho0 * US%R_to_kg_m3*US%L_T_to_m_s**2*US%Z_to_L*tau_mag)
endif
enddo ; enddo
else ! C-grid wind stresses.
Expand All @@ -1062,11 +1064,11 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
gustiness = CS%gust_const
if (CS%read_gust_2d) gustiness = CS%gust(i,j)

if (do_ustar) ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * tau_mag)
if (do_ustar) ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * US%R_to_kg_m3*US%L_T_to_m_s**2*US%Z_to_L*tau_mag)
if (CS%answers_2018) then
if (do_gustless) gustless_ustar(i,j) = US%m_to_Z*US%T_to_s * sqrt(tau_mag / CS%Rho0)
if (do_gustless) gustless_ustar(i,j) = sqrt(US%R_to_kg_m3*US%L_to_Z*tau_mag / CS%Rho0)
else
if (do_gustless) gustless_ustar(i,j) = sqrt(Irho0 * tau_mag)
if (do_gustless) gustless_ustar(i,j) = sqrt(Irho0 * US%R_to_kg_m3*US%L_T_to_m_s**2*US%Z_to_L*tau_mag)
endif
enddo ; enddo
endif ! endif for wind friction velocity fields
Expand Down Expand Up @@ -1132,12 +1134,15 @@ subroutine apply_force_adjustments(G, CS, Time, forces)
type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces

! Local variables
real, dimension(SZI_(G),SZJ_(G)) :: tempx_at_h ! Delta to zonal wind stress at h points [Pa]
real, dimension(SZI_(G),SZJ_(G)) :: tempy_at_h ! Delta to meridional wind stress at h points [Pa]
real, dimension(SZI_(G),SZJ_(G)) :: tempx_at_h ! Delta to zonal wind stress at h points [R Z L T-2 ~> Pa]
real, dimension(SZI_(G),SZJ_(G)) :: tempy_at_h ! Delta to meridional wind stress at h points [R Z L T-2 ~> Pa]

integer :: isc, iec, jsc, jec, i, j
real :: dLonDx, dLonDy, rDlon, cosA, sinA, zonal_tau, merid_tau
logical :: overrode_x, overrode_y
type(unit_scale_type), pointer :: US => NULL() !< A dimensional unit scaling type

US => G%US

isc = G%isc; iec = G%iec ; jsc = G%jsc; jec = G%jec

Expand All @@ -1160,8 +1165,8 @@ subroutine apply_force_adjustments(G, CS, Time, forces)
if (rDlon > 0.) rDlon = 1. / rDlon
cosA = dLonDx * rDlon
sinA = dLonDy * rDlon
zonal_tau = tempx_at_h(i,j)
merid_tau = tempy_at_h(i,j)
zonal_tau = US%kg_m3_to_R*US%m_s_to_L_T**2*US%L_to_Z * tempx_at_h(i,j)
merid_tau = US%kg_m3_to_R*US%m_s_to_L_T**2*US%L_to_Z * tempy_at_h(i,j)
tempx_at_h(i,j) = cosA * zonal_tau - sinA * merid_tau
tempy_at_h(i,j) = sinA * zonal_tau + cosA * merid_tau
enddo ; enddo
Expand Down
30 changes: 17 additions & 13 deletions config_src/ice_solo_driver/MOM_surface_forcing.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -389,8 +389,8 @@ subroutine wind_forcing_2gyre(sfc_state, forces, day, G, CS)
PI = 4.0*atan(1.0)

do j=js,je ; do I=Isq,Ieq
forces%taux(I,j) = 0.1*(1.0 - cos(2.0*PI*(G%geoLatCu(I,j)-CS%South_lat) / &
CS%len_lat))
forces%taux(I,j) = 0.1*US%kg_m3_to_R*US%m_s_to_L_T**2*US%L_to_Z * &
(1.0 - cos(2.0*PI*(G%geoLatCu(I,j)-CS%South_lat) / CS%len_lat))
enddo ; enddo

do J=Jsq,Jeq ; do i=is,ie
Expand Down Expand Up @@ -426,7 +426,8 @@ subroutine wind_forcing_1gyre(sfc_state, forces, day, G, CS)
PI = 4.0*atan(1.0)

do j=js,je ; do I=Isq,Ieq
forces%taux(I,j) =-0.2*cos(PI*(G%geoLatCu(I,j)-CS%South_lat)/CS%len_lat)
forces%taux(I,j) = -0.2*US%kg_m3_to_R*US%m_s_to_L_T**2*US%L_to_Z * &
cos(PI*(G%geoLatCu(I,j)-CS%South_lat)/CS%len_lat)
enddo ; enddo

do J=Jsq,Jeq ; do i=is,ie
Expand Down Expand Up @@ -464,9 +465,9 @@ subroutine wind_forcing_gyres(sfc_state, forces, day, G, US, CS)

do j=jsd,jed ; do I=IsdB,IedB
y = (G%geoLatCu(I,j)-CS%South_lat)/CS%len_lat
forces%taux(I,j) = CS%gyres_taux_const + &
forces%taux(I,j) = US%kg_m3_to_R*US%m_s_to_L_T**2*US%L_to_Z * (CS%gyres_taux_const + &
( CS%gyres_taux_sin_amp*sin(CS%gyres_taux_n_pis*PI*y) &
+ CS%gyres_taux_cos_amp*cos(CS%gyres_taux_n_pis*PI*y) )
+ CS%gyres_taux_cos_amp*cos(CS%gyres_taux_n_pis*PI*y) ))
enddo ; enddo

do J=JsdB,JedB ; do i=isd,ied
Expand All @@ -477,7 +478,7 @@ subroutine wind_forcing_gyres(sfc_state, forces, day, G, US, CS)
do j=js,je ; do i=is,ie
forces%ustar(i,j) = US%m_to_Z*US%T_to_s * sqrt(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 + (CS%gust_const/CS%Rho0))
forces%taux(i,j)*forces%taux(i,j)))* US%R_to_kg_m3*US%L_T_to_m_s**2*US%Z_to_L /CS%Rho0 + (CS%gust_const/CS%Rho0))
enddo ; enddo

call callTree_leave("wind_forcing_gyres")
Expand Down Expand Up @@ -528,10 +529,12 @@ subroutine wind_forcing_from_file(sfc_state, forces, day, G, US, CS)

call pass_vector(temp_x, temp_y, G%Domain, To_All, AGRID)
do j=js,je ; do I=Isq,Ieq
forces%taux(I,j) = 0.5 * CS%wind_scale * (temp_x(i,j) + temp_x(i+1,j))
forces%taux(I,j) = 0.5 * CS%wind_scale * US%kg_m3_to_R*US%m_s_to_L_T**2*US%L_to_Z * &
(temp_x(i,j) + temp_x(i+1,j))
enddo ; enddo
do J=Jsq,Jeq ; do i=is,ie
forces%tauy(i,J) = 0.5 * CS%wind_scale * (temp_y(i,j) + temp_y(i,j+1))
forces%tauy(i,J) = 0.5 * CS%wind_scale * US%kg_m3_to_R*US%m_s_to_L_T**2*US%L_to_Z * &
(temp_y(i,j) + temp_y(i,j+1))
enddo ; enddo

if (CS%read_gust_2d) then
Expand All @@ -548,7 +551,8 @@ subroutine wind_forcing_from_file(sfc_state, forces, day, G, US, CS)
case ("C")
call MOM_read_vector(filename,CS%stress_x_var, CS%stress_y_var, &
forces%taux(:,:), forces%tauy(:,:), &
G%Domain, timelevel=time_lev)
G%Domain, timelevel=time_lev, &
scale=US%kg_m3_to_R*US%m_s_to_L_T**2*US%L_to_Z)
if (CS%wind_scale /= 1.0) then
do j=js,je ; do I=Isq,Ieq
forces%taux(I,j) = CS%wind_scale * forces%taux(I,j)
Expand All @@ -561,15 +565,15 @@ subroutine wind_forcing_from_file(sfc_state, forces, day, G, US, CS)
call pass_vector(forces%taux, forces%tauy, G%Domain, To_All)
if (CS%read_gust_2d) then
do j=js, je ; do i=is, ie
forces%ustar(i,j) = US%m_to_Z*US%T_to_s * sqrt((sqrt(0.5*((forces%tauy(i,j-1)**2 + &
forces%tauy(i,j)**2) + (forces%taux(i-1,j)**2 + &
forces%ustar(i,j) = US%m_to_Z*US%T_to_s * sqrt((US%R_to_kg_m3*US%L_T_to_m_s**2*US%Z_to_L * &
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%gust(i,j)) / CS%Rho0 )
enddo ; enddo
else
do j=js, je ; do i=is, ie
forces%ustar(i,j) = US%m_to_Z*US%T_to_s * sqrt(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 + (CS%gust_const/CS%Rho0))
forces%tauy(i,j)**2) + (forces%taux(i-1,j)**2 + forces%taux(i,j)**2))) * &
US%R_to_kg_m3*US%L_T_to_m_s**2*US%Z_to_L/CS%Rho0 + (CS%gust_const/CS%Rho0))
enddo ; enddo
endif
case default
Expand Down
9 changes: 5 additions & 4 deletions config_src/ice_solo_driver/user_surface_forcing.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -91,7 +91,7 @@ module user_surface_forcing

contains

!> This subroutine sets the surface wind stresses, forces%taux and forces%tauy, in [Pa].
!> This subroutine sets the surface wind stresses, forces%taux and forces%tauy, in [R Z L T-2 ~> Pa].
!! These are the stresses in the direction of the model grid (i.e. the same
!! direction as the u- and v- velocities).
subroutine USER_wind_forcing(sfc_state, forces, day, G, US, CS)
Expand All @@ -104,7 +104,7 @@ subroutine USER_wind_forcing(sfc_state, forces, day, G, US, CS)
type(user_surface_forcing_CS), pointer :: CS !< A pointer to the control structure returned
!! by a previous call to user_surface_forcing_init

! This subroutine sets the surface wind stresses, forces%taux and forces%tauy [Pa].
! This subroutine sets the surface wind stresses, forces%taux and forces%tauy [R Z L T-2 ~> Pa].
! In addition, this subroutine can be used to set the surface friction velocity,
! forces%ustar [Z T-1 ~> m s-1], which is needed with a bulk mixed layer.

Expand All @@ -130,7 +130,8 @@ subroutine USER_wind_forcing(sfc_state, forces, day, G, US, CS)
! The i-loop extends to is-1 so that taux can be used later in the
! calculation of ustar - otherwise the lower bound would be Isq.
do j=js,je ; do I=is-1,Ieq
forces%taux(I,j) = G%mask2dCu(I,j) * 0.0 ! Change this to the desired expression.
! Change this to the desired expression.
forces%taux(I,j) = G%mask2dCu(I,j) * 0.0*US%kg_m3_to_R*US%m_s_to_L_T**2*US%L_to_Z
enddo ; enddo
do J=js-1,Jeq ; do i=is,ie
forces%tauy(i,J) = G%mask2dCv(i,J) * 0.0 ! Change this to the desired expression.
Expand All @@ -140,7 +141,7 @@ 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%ustar(i,j) = US%m_to_Z*US%T_to_s * G%mask2dT(i,j) * sqrt(CS%gust_const/CS%Rho0 + &
sqrt(0.5*(forces%taux(I-1,j)**2 + forces%taux(I,j)**2) + &
US%R_to_kg_m3*US%L_T_to_m_s**2*US%Z_to_L*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

Expand Down
Loading

0 comments on commit 1cbf498

Please sign in to comment.