Skip to content

Commit

Permalink
Rescaled internal MOM_CVMix variables
Browse files Browse the repository at this point in the history 
  Applied dimensional rescaling to many of the internal calculations in the 4
MOM_CVMix files, although calls to external CVMix routines still use the
original MKS units.  These changes include rescaling of the input and output
variables associated with the calculate_density routines.  All answers are
bitwise identical.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Apr 17, 2020
1 parent 4f42d72 commit 9b53927
Show file tree
Hide file tree
Showing 4 changed files with 75 additions and 142 deletions.
110 changes: 16 additions & 94 deletions src/parameterizations/vertical/MOM_CVMix_KPP.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -148,7 +148,7 @@ module MOM_CVMix_KPP
real, allocatable, dimension(:,:) :: kOBL !< Level (+fraction) of OBL extent
real, allocatable, dimension(:,:) :: OBLdepthprev !< previous Depth (positive) of OBL [m]
real, allocatable, dimension(:,:) :: La_SL !< Langmuir number used in KPP
real, allocatable, dimension(:,:,:) :: dRho !< Bulk difference in density [kg m-3]
real, allocatable, dimension(:,:,:) :: dRho !< Bulk difference in density [R ~> kg m-3]
real, allocatable, dimension(:,:,:) :: Uz2 !< Square of bulk difference in resolved velocity [m2 s-2]
real, allocatable, dimension(:,:,:) :: BulkRi !< Bulk Richardson number for each layer (dimensionless)
real, allocatable, dimension(:,:,:) :: sigma !< Sigma coordinate (dimensionless)
Expand Down Expand Up @@ -188,7 +188,7 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves)
type(wave_parameters_CS), optional, pointer :: Waves !< Wave CS

! Local variables
#include "version_variable.h"
# include "version_variable.h"
character(len=40) :: mdl = 'MOM_CVMix_KPP' !< name of this module
character(len=20) :: string !< local temporary string
logical :: CS_IS_ONE=.false. !< Logical for setting Cs based on Non-local
Expand Down Expand Up @@ -475,7 +475,8 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves)
cmor_units='m', cmor_standard_name='Ocean Mixed Layer Thickness Defined by Mixing Scheme')
endif
CS%id_BulkDrho = register_diag_field('ocean_model', 'KPP_BulkDrho', diag%axesTL, Time, &
'Bulk difference in density used in Bulk Richardson number, as used by [CVMix] KPP', 'kg/m3')
'Bulk difference in density used in Bulk Richardson number, as used by [CVMix] KPP', &
'kg/m3', conversion=US%R_to_kg_m3)
CS%id_BulkUz2 = register_diag_field('ocean_model', 'KPP_BulkUz2', diag%axesTL, Time, &
'Square of bulk difference in resolved velocity used in Bulk Richardson number via [CVMix] KPP', 'm2/s2')
CS%id_BulkRi = register_diag_field('ocean_model', 'KPP_BulkRi', diag%axesTL, Time, &
Expand Down Expand Up @@ -908,20 +909,21 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, EOS, uStar, buoyF
real, dimension( G%ke+1 ) :: iFaceHeight ! Interface heights referenced to surface [m] (negative in ocean)
real, dimension( G%ke+1 ) :: N2_1d ! Brunt-Vaisala frequency squared, at interfaces [s-2]
real, dimension( G%ke ) :: Ws_1d ! Profile of vertical velocity scale for scalars [m s-1]
real, dimension( G%ke ) :: deltaRho ! delta Rho in numerator of Bulk Ri number
real, dimension( G%ke ) :: deltaRho ! delta Rho in numerator of Bulk Ri number [R ~> kg m-3]
real, dimension( G%ke ) :: deltaU2 ! square of delta U (shear) in denominator of Bulk Ri [m2 s-2]
real, dimension( G%ke ) :: surfBuoyFlux2
real, dimension( G%ke ) :: BulkRi_1d ! Bulk Richardson number for each layer

! for EOS calculation
real, dimension( 3*G%ke ) :: rho_1D
real, dimension( 3*G%ke ) :: pres_1D
real, dimension( 3*G%ke ) :: rho_1D ! A column of densities [R ~> kg m-3]
real, dimension( 3*G%ke ) :: pres_1D ! A column of pressures [R L2 T-2 ~> Pa]
real, dimension( 3*G%ke ) :: Temp_1D
real, dimension( 3*G%ke ) :: Salt_1D

real :: surfFricVel, surfBuoyFlux, Coriolis
real :: GoRho ! Gravitational acceleration divided by density in MKS units [m4 s-2]
real :: pRef, rho1, rhoK, Uk, Vk, sigma, sigmaRatio
real :: GoRho ! Gravitational acceleration divided by density in MKS units [m R-1 s-2 ~> m4 kg-1 s-2]
real :: pRef ! The interface pressure [R L2 T-2 ~> Pa]
real :: rho1, rhoK, Uk, Vk, sigma, sigmaRatio

real :: zBottomMinusOffset ! Height of bottom plus a little bit [m]
real :: SLdepth_0d ! Surface layer depth = surf_layer_ext*OBLdepth.
Expand Down Expand Up @@ -958,7 +960,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, EOS, uStar, buoyF
#endif

! some constants
GoRho = GV%mks_g_Earth / (US%R_to_kg_m3*GV%Rho0)
GoRho = GV%mks_g_Earth / GV%Rho0
buoy_scale = US%L_to_m**2*US%s_to_T**3

! loop over horizontal points on processor
Expand Down Expand Up @@ -1084,9 +1086,9 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, EOS, uStar, buoyF
Salt_1D(kk+2) = Salt(i,j,k)
Salt_1D(kk+3) = Salt(i,j,km1)

! pRef is pressure at interface between k and km1.
! pRef is pressure at interface between k and km1 [R L2 T-2 ~> Pa].
! iterate pRef for next pass through k-loop.
pRef = pRef + GV%H_to_Pa * h(i,j,k)
pRef = pRef + (GV%g_Earth * GV%H_to_RZ) * h(i,j,k)

! this difference accounts for penetrating SW
surfBuoyFlux2(k) = buoy_scale * (buoyFlux(i,j,1) - buoyFlux(i,j,k+1))
Expand All @@ -1102,7 +1104,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, EOS, uStar, buoyF


! compute in-situ density
call calculate_density(Temp_1D, Salt_1D, pres_1D, rho_1D, 1, 3*G%ke, EOS)
call calculate_density(Temp_1D, Salt_1D, pres_1D, rho_1D, EOS, US)

! N2 (can be negative) and N (non-negative) on interfaces.
! deltaRho is non-local rho difference used for bulk Richardson number.
Expand Down Expand Up @@ -1215,86 +1217,6 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, EOS, uStar, buoyF
CS%OBLdepth(i,j) = min( CS%OBLdepth(i,j), -iFaceHeight(G%ke+1) ) ! no deeper than bottom
CS%kOBL(i,j) = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) )

!*************************************************************************
! smg: remove code below

! Following "correction" step has been found to be unnecessary.
! Code should be removed after further testing.
! BGR: 03/15/2018-> Restructured code (Vt2 changed to compute from call in MOM_CVMix_KPP now)
! I have not taken this restructuring into account here.
! Do we ever run with correctSurfLayerAvg?
! smg's suggested testing and removal is advised, in the meantime
! I have added warning if correctSurfLayerAvg is attempted.
! if (CS%correctSurfLayerAvg) then

! SLdepth_0d = CS%surf_layer_ext * CS%OBLdepth(i,j)
! hTot = h(i,j,1)
! surfTemp = Temp(i,j,1) ; surfHtemp = surfTemp * hTot
! surfSalt = Salt(i,j,1) ; surfHsalt = surfSalt * hTot
! surfU = 0.5*US%L_T_to_m_s*(u(i,j,1)+u(i-1,j,1)) ; surfHu = surfU * hTot
! surfV = 0.5*US%L_T_to_m_s*(v(i,j,1)+v(i,j-1,1)) ; surfHv = surfV * hTot
! pRef = 0.0

! do k = 2, G%ke

! ! Recalculate differences with surface layer
! Uk = 0.5*US%L_T_to_m_s*(u(i,j,k)+u(i-1,j,k)) - surfU
! Vk = 0.5*US%L_T_to_m_s*(v(i,j,k)+v(i,j-1,k)) - surfV
! deltaU2(k) = Uk**2 + Vk**2
! pRef = pRef + GV%H_to_Pa * h(i,j,k)
! call calculate_density(surfTemp, surfSalt, pRef, rho1, EOS)
! call calculate_density(Temp(i,j,k), Salt(i,j,k), pRef, rhoK, EOS)
! deltaRho(k) = rhoK - rho1

! ! Surface layer averaging (needed for next k+1 iteration of this loop)
! if (hTot < SLdepth_0d) then
! delH = min( max(0., SLdepth_0d - hTot), h(i,j,k)*GV%H_to_m )
! hTot = hTot + delH
! surfHtemp = surfHtemp + Temp(i,j,k) * delH ; surfTemp = surfHtemp / hTot
! surfHsalt = surfHsalt + Salt(i,j,k) * delH ; surfSalt = surfHsalt / hTot
! surfHu = surfHu + 0.5*US%L_T_to_m_s*(u(i,j,k)+u(i-1,j,k)) * delH ; surfU = surfHu / hTot
! surfHv = surfHv + 0.5*US%L_T_to_m_s*(v(i,j,k)+v(i,j-1,k)) * delH ; surfV = surfHv / hTot
! endif

! enddo

! BulkRi_1d = CVMix_kpp_compute_bulk_Richardson( &
! cellHeight(1:G%ke), & ! Depth of cell center [m]
! GoRho*deltaRho, & ! Bulk buoyancy difference, Br-B(z) [s-1]
! deltaU2, & ! Square of resolved velocity difference [m2 s-2]
! ws_cntr=Ws_1d, & ! Turbulent velocity scale profile [m s-1]
! N_iface=CS%N ) ! Buoyancy frequency [s-1]

! surfBuoyFlux = buoy_scale*buoyFlux(i,j,1) ! This is only used in kpp_compute_OBL_depth to limit
! ! h to Monin-Obukov (default is false, ie. not used)

! call CVMix_kpp_compute_OBL_depth( &
! BulkRi_1d, & ! (in) Bulk Richardson number
! iFaceHeight, & ! (in) Height of interfaces [m]
! CS%OBLdepth(i,j), & ! (out) OBL depth [m]
! CS%kOBL(i,j), & ! (out) level (+fraction) of OBL extent
! zt_cntr=cellHeight, & ! (in) Height of cell centers [m]
! surf_fric=surfFricVel, & ! (in) Turbulent friction velocity at surface [m s-1]
! surf_buoy=surfBuoyFlux, & ! (in) Buoyancy flux at surface [m2 s-3]
! Coriolis=Coriolis, & ! (in) Coriolis parameter [s-1]
! CVMix_kpp_params_user=CS%KPP_params ) ! KPP parameters

! if (CS%deepOBLoffset>0.) then
! zBottomMinusOffset = iFaceHeight(G%ke+1) + min(CS%deepOBLoffset,-0.1*iFaceHeight(G%ke+1))
! CS%OBLdepth(i,j) = min( CS%OBLdepth(i,j), -zBottomMinusOffset )
! CS%kOBL(i,j) = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) )
! endif

! ! apply some constraints on OBLdepth
! if(CS%fixedOBLdepth) CS%OBLdepth(i,j) = CS%fixedOBLdepth_value
! CS%OBLdepth(i,j) = max( CS%OBLdepth(i,j), -iFaceHeight(2) ) ! no shallower than top layer
! CS%OBLdepth(i,j) = min( CS%OBLdepth(i,j), -iFaceHeight(G%ke+1) ) ! no deep than bottom
! CS%kOBL(i,j) = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) )

! endif ! endif for "correction" step

! smg: remove code above
! **********************************************************************

! recompute wscale for diagnostics, now that we in fact know boundary layer depth
!BGR consider if LTEnhancement is wanted for diagnostics
Expand Down Expand Up @@ -1359,7 +1281,7 @@ subroutine KPP_smooth_BLD(CS,G,GV,h)
real :: wc, ww, we, wn, ws ! averaging weights for smoothing
real :: dh ! The local thickness used for calculating interface positions [m]
real :: hcorr ! A cumulative correction arising from inflation of vanished layers [m]
real :: pref
!### real :: pref
integer :: i, j, k, s

do s=1,CS%n_smooth
Expand All @@ -1378,7 +1300,7 @@ subroutine KPP_smooth_BLD(CS,G,GV,h)
if (G%mask2dT(i,j)==0.) cycle

iFaceHeight(1) = 0.0 ! BBL is all relative to the surface
pRef = 0.
!### pRef = 0.
hcorr = 0.
do k=1,G%ke

Expand Down
19 changes: 11 additions & 8 deletions src/parameterizations/vertical/MOM_CVMix_conv.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -168,11 +168,14 @@ subroutine calculate_CVMix_conv(h, tv, G, GV, US, CS, hbl)
real, dimension(SZK_(G)+1) :: iFaceHeight !< Height of interfaces [m]
real, dimension(SZK_(G)) :: cellHeight !< Height of cell centers [m]
integer :: kOBL !< level of OBL extent
real :: g_o_rho0 ! Gravitational acceleration divided by density in MKS units [m4 s-2]
real :: pref, rhok, rhokm1, dz, dh, hcorr
real :: g_o_rho0 ! Gravitational acceleration divided by density times unit convserion factors
! [Z s-2 R-1 ~> m4 s-2 kg-1]
real :: pref ! Interface pressures [R L2 T-2 ~> Pa]
real :: rhok, rhokm1 ! In situ densities of the layers above and below at the interface pressure [R ~> kg m-3]
real :: dz, dh, hcorr
integer :: i, j, k

g_o_rho0 = GV%mks_g_Earth / (US%R_to_kg_m3*GV%Rho0)
g_o_rho0 = US%L_to_Z**2*US%s_to_T**2 * GV%g_Earth / GV%Rho0

! initialize dummy variables
rho_lwr(:) = 0.0; rho_1d(:) = 0.0
Expand All @@ -196,12 +199,12 @@ subroutine calculate_CVMix_conv(h, tv, G, GV, US, CS, hbl)
! Compute Brunt-Vaisala frequency (static stability) on interfaces
do k=2,G%ke

! pRef is pressure at interface between k and km1.
pRef = pRef + GV%H_to_Pa * h(i,j,k)
call calculate_density(tv%t(i,j,k), tv%s(i,j,k), pref, rhok, tv%eqn_of_state)
call calculate_density(tv%t(i,j,k-1), tv%s(i,j,k-1), pref, rhokm1, tv%eqn_of_state)
! pRef is pressure at interface between k and km1 [R L2 T-2 ~> Pa].
pRef = pRef + (GV%H_to_RZ*GV%g_Earth) * h(i,j,k)
call calculate_density(tv%t(i,j,k), tv%s(i,j,k), pRef, rhok, tv%eqn_of_state, US=US)
call calculate_density(tv%t(i,j,k-1), tv%s(i,j,k-1), pRef, rhokm1, tv%eqn_of_state, US=US)

dz = ((0.5*(h(i,j,k-1) + h(i,j,k))+GV%H_subroundoff)*GV%H_to_m)
dz = ((0.5*(h(i,j,k-1) + h(i,j,k))+GV%H_subroundoff)*GV%H_to_Z)
CS%N2(i,j,k) = g_o_rho0 * (rhok - rhokm1) / dz ! Can be negative

enddo
Expand Down
34 changes: 16 additions & 18 deletions src/parameterizations/vertical/MOM_CVMix_ddiff.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -182,13 +182,13 @@ subroutine compute_ddiff_coeffs(h, tv, G, GV, US, j, Kd_T, Kd_S, CS)
! Local variables
real, dimension(SZK_(G)) :: &
cellHeight, & !< Height of cell centers [m]
dRho_dT, & !< partial derivatives of density wrt temp [kg m-3 degC-1]
dRho_dS, & !< partial derivatives of density wrt saln [kg m-3 ppt-1]
pres_int, & !< pressure at each interface [Pa]
dRho_dT, & !< partial derivatives of density wrt temp [R degC-1 ~> kg m-3 degC-1]
dRho_dS, & !< partial derivatives of density wrt saln [R ppt-1 ~> kg m-3 ppt-1]
pres_int, & !< pressure at each interface [R L2 T-2 ~> Pa]
temp_int, & !< temp and at interfaces [degC]
salt_int, & !< salt at at interfaces [ppt]
alpha_dT, & !< alpha*dT across interfaces
beta_dS, & !< beta*dS across interfaces
alpha_dT, & !< alpha*dT across interfaces [kg m-3]
beta_dS, & !< beta*dS across interfaces [kg m-3]
dT, & !< temp. difference between adjacent layers [degC]
dS !< salt difference between adjacent layers [ppt]
real, dimension(SZK_(G)+1) :: &
Expand All @@ -197,7 +197,7 @@ subroutine compute_ddiff_coeffs(h, tv, G, GV, US, j, Kd_T, Kd_S, CS)

real, dimension(SZK_(G)+1) :: iFaceHeight !< Height of interfaces [m]
integer :: kOBL !< level of OBL extent
real :: pref, g_o_rho0, rhok, rhokm1, dz, dh, hcorr
real :: dh, hcorr
integer :: i, k

! initialize dummy variables
Expand All @@ -219,31 +219,29 @@ subroutine compute_ddiff_coeffs(h, tv, G, GV, US, j, Kd_T, Kd_S, CS)
! skip calling at land points
if (G%mask2dT(i,j) == 0.) cycle

pRef = 0.
pres_int(1) = pRef
pres_int(1) = 0.
! we don't have SST and SSS, so let's use values at top-most layer
temp_int(1) = TV%T(i,j,1); salt_int(1) = TV%S(i,j,1)
do k=2,G%ke
do K=2,G%ke
! pressure at interface
pres_int(k) = pRef + GV%H_to_Pa * h(i,j,k-1)
pres_int(K) = pres_int(K-1) + (GV%g_Earth * GV%H_to_RZ) * h(i,j,k-1)
! temp and salt at interface
! for temp: (t1*h1 + t2*h2)/(h1+h2)
temp_int(k) = (TV%T(i,j,k-1)*h(i,j,k-1) + TV%T(i,j,k)*h(i,j,k))/(h(i,j,k-1)+h(i,j,k))
salt_int(k) = (TV%S(i,j,k-1)*h(i,j,k-1) + TV%S(i,j,k)*h(i,j,k))/(h(i,j,k-1)+h(i,j,k))
temp_int(K) = (TV%T(i,j,k-1)*h(i,j,k-1) + TV%T(i,j,k)*h(i,j,k))/(h(i,j,k-1)+h(i,j,k))
salt_int(K) = (TV%S(i,j,k-1)*h(i,j,k-1) + TV%S(i,j,k)*h(i,j,k))/(h(i,j,k-1)+h(i,j,k))
! dT and dS
dT(k) = (TV%T(i,j,k-1)-TV%T(i,j,k))
dS(k) = (TV%S(i,j,k-1)-TV%S(i,j,k))
pRef = pRef + GV%H_to_Pa * h(i,j,k-1)
dT(K) = (TV%T(i,j,k-1)-TV%T(i,j,k))
dS(K) = (TV%S(i,j,k-1)-TV%S(i,j,k))
enddo ! k-loop finishes

call calculate_density_derivs(temp_int, salt_int, pres_int, drho_dT, drho_dS, tv%eqn_of_state)
call calculate_density_derivs(temp_int, salt_int, pres_int, drho_dT, drho_dS, tv%eqn_of_state, US)

! The "-1.0" below is needed so that the following criteria is satisfied:
! if ((alpha_dT > beta_dS) .and. (beta_dS > 0.0)) then "salt finger"
! if ((alpha_dT < 0.) .and. (beta_dS < 0.) .and. (alpha_dT > beta_dS)) then "diffusive convection"
do k=1,G%ke
alpha_dT(k) = -1.0*drho_dT(k) * dT(k)
beta_dS(k) = drho_dS(k) * dS(k)
alpha_dT(k) = -1.0*US%R_to_kg_m3*drho_dT(k) * dT(k)
beta_dS(k) = US%R_to_kg_m3*drho_dS(k) * dS(k)
enddo

if (CS%id_R_rho > 0.0) then
Expand Down
Loading

0 comments on commit 9b53927

Please sign in to comment.