Skip to content

Commit

Permalink
(*)Avoid negative thicknesses in mixed_layer_restrat
Browse files Browse the repository at this point in the history 
  Enforce a minimum thickness of 0.5*Angstrom in the mixed_layer_restrat
routines.  The streamfunctions in these routines attempt to limit the
thicknesses to exceed Angstrom, but they can be less than this due to roundoff.
The new limit prevents thicknesses from becoming negative when Angstrom is set
to 0, but should not change any answers for test cases with positive values of
Angstrom.  Also added some comments describing variables and their units and
simplified the OMP directives.

  Also corrected error messages in MOM_diabatic_aux.F90 to identify the file or
module where these messages come from, and modified an error message in
applyTracerBoundaryFluxesInOut so that it is written if the localized fault does
not happen to occur on the root PE.

  All answers in the existing MOM6-examples regression suite are bitwise
identical.
  • Loading branch information
@Hallberg-NOAA @marshallward
Hallberg-NOAA authored and marshallward committed Mar 7, 2022
1 parent 9a01cd5 commit 4640461
Show file tree
Hide file tree
Showing 3 changed files with 50 additions and 49 deletions.
85 changes: 43 additions & 42 deletions src/parameterizations/lateral/MOM_mixed_layer_restrat.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -156,37 +156,42 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var
real :: u_star ! surface friction velocity, interpolated to velocity points [Z T-1 ~> m s-1].
real :: mom_mixrate ! rate at which momentum is homogenized within mixed layer [T-1 ~> s-1]
real :: timescale ! mixing growth timescale [T ~> s]
real :: h_min ! The minimum layer thickness [H ~> m or kg m-2]. h_min could be 0.
real :: h_neglect ! tiny thickness usually lost in roundoff so can be neglected [H ~> m or kg m-2]
real :: dz_neglect ! A tiny thickness that is usually lost in roundoff so can be neglected [Z ~> m]
real :: I4dt ! 1/(4 dt) [T-1 ~> s-1]
real :: Ihtot,Ihtot_slow! Inverses of the total mixed layer thickness [H-1 ~> m-1 or m2 kg-1]
real :: a(SZK_(GV)) ! A non-dimensional value relating the overall flux
! magnitudes (uDml & vDml) to the realized flux in a
! layer. The vertical sum of a() through the pieces of
! layer [nondim]. The vertical sum of a() through the pieces of
! the mixed layer must be 0.
real :: b(SZK_(GV)) ! As for a(k) but for the slow-filtered MLD
real :: b(SZK_(GV)) ! As for a(k) but for the slow-filtered MLD [nondim]
real :: uDml(SZIB_(G)) ! The zonal and meridional volume fluxes in the upper
real :: vDml(SZI_(G)) ! half of the mixed layer [H L2 T-1 ~> m3 s-1 or kg s-1].
real :: uDml_slow(SZIB_(G)) ! The zonal and meridional volume fluxes in the upper
real :: vDml_slow(SZI_(G)) ! half of the mixed layer [H L2 T-1 ~> m3 s-1 or kg s-1].
real :: utimescale_diag(SZIB_(G),SZJ_(G)) ! restratification timescales in the zonal and
real :: vtimescale_diag(SZI_(G),SZJB_(G)) ! meridional directions [T ~> s], stored in 2-D arrays
! for diagnostic purposes.
real :: uDml_diag(SZIB_(G),SZJ_(G)), vDml_diag(SZI_(G),SZJB_(G))
real :: uDml_diag(SZIB_(G),SZJ_(G)) ! A 2D copy of uDml for diagnostics [H L2 T-1 ~> m3 s-1 or kg s-1]
real :: vDml_diag(SZI_(G),SZJB_(G)) ! A 2D copy of vDml for diagnostics [H L2 T-1 ~> m3 s-1 or kg s-1]
real, dimension(SZI_(G)) :: rhoSurf, deltaRhoAtKm1, deltaRhoAtK ! Densities [R ~> kg m-3]
real, dimension(SZI_(G)) :: dK, dKm1 ! Depths of layer centers [H ~> m or kg m-2].
real, dimension(SZI_(G)) :: pRef_MLD ! A reference pressure for calculating the mixed layer
! densities [R L2 T-2 ~> Pa].
real, dimension(SZI_(G)) :: rhoAtK, rho1, d1, pRef_N2 ! Used for N2
real :: aFac, bFac ! Nondimensional ratios [nondim]
real :: ddRho ! A density difference [R ~> kg m-3]
real :: hAtVel, zpa, zpb, dh, res_scaling_fac
real :: aFac, bFac ! Nondimensional ratios [nondim]
real :: ddRho ! A density difference [R ~> kg m-3]
real :: hAtVel ! Thickness at the velocity points [H ~> m or kg m-2]
real :: zpa ! Fractional position within the mixed layer of the interface above a layer [nondim]
real :: zpb ! Fractional position within the mixed layer of the interface below a layer [nondim]
real :: dh ! Portion of the layer thickness that is in the mixed layer [H ~> m or kg m-2]
real :: res_scaling_fac ! The resolution-dependent scaling factor [nondim]
real :: I_LFront ! The inverse of the frontal length scale [L-1 ~> m-1]
logical :: line_is_empty, keep_going, res_upscale
integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz

real :: PSI, PSI1, z, BOTTOP, XP, DD ! For the following statement functions
real :: PSI, PSI1, z, BOTTOP, XP, DD ! For the following statement functions [nondim]
! Stream function as a function of non-dimensional position within mixed-layer (F77 statement function)
!PSI1(z) = max(0., (1. - (2.*z+1.)**2 ) )
PSI1(z) = max(0., (1. - (2.*z+1.)**2 ) * (1. + (5./21.)*(2.*z+1.)**2) )
Expand All @@ -198,6 +203,8 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB

h_min = 0.5*GV%Angstrom_H ! This should be GV%Angstrom_H, but that value would change answers.

if (.not.associated(tv%eqn_of_state)) call MOM_error(FATAL, "MOM_mixedlayer_restrat: "// &
"An equation of state must be used with this module.")
if (.not. allocated(VarMix%Rd_dx_h) .and. CS%front_length > 0.) &
Expand Down Expand Up @@ -299,16 +306,11 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var

p0(:) = 0.0
EOSdom(:) = EOS_domain(G%HI, halo=1)
!$OMP parallel default(none) shared(is,ie,js,je,G,GV,US,htot_fast,Rml_av_fast,tv,p0,h,h_avail,&
!$OMP h_neglect,g_Rho0,I4dt,CS,uhml,uhtr,dt,vhml,vhtr,EOSdom, &
!$OMP utimescale_diag,vtimescale_diag,forces,dz_neglect, &
!$OMP htot_slow,MLD_slow,Rml_av_slow,VarMix,I_LFront, &
!$OMP res_upscale, nz,MLD_fast,uDml_diag,vDml_diag) &
!$OMP private(rho_ml,h_vel,u_star,absf,mom_mixrate,timescale, &
!$OMP line_is_empty, keep_going,res_scaling_fac, &
!$OMP a,IhTot,b,Ihtot_slow,zpb,hAtVel,zpa,dh) &
!$OMP firstprivate(uDml,vDml,uDml_slow,vDml_slow)
!$OMP do
!$OMP parallel default(shared) private(rho_ml,h_vel,u_star,absf,mom_mixrate,timescale, &
!$OMP line_is_empty, keep_going,res_scaling_fac, &
!$OMP a,IhTot,b,Ihtot_slow,zpb,hAtVel,zpa,dh) &
!$OMP firstprivate(uDml,vDml,uDml_slow,vDml_slow)
!$OMP do
do j=js-1,je+1
do i=is-1,ie+1
htot_fast(i,j) = 0.0 ; Rml_av_fast(i,j) = 0.0
Expand Down Expand Up @@ -359,7 +361,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var
! 2. Add exponential tail to stream-function?

! U - Component
!$OMP do
!$OMP do
do j=js,je ; do I=is-1,ie
u_star = 0.5*(forces%ustar(i,j) + forces%ustar(i+1,j))
absf = 0.5*(abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I,J)))
Expand Down Expand Up @@ -435,7 +437,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var
enddo ; enddo

! V- component
!$OMP do
!$OMP do
do J=js-1,je ; do i=is,ie
u_star = 0.5*(forces%ustar(i,j) + forces%ustar(i,j+1))
absf = 0.5*(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J)))
Expand Down Expand Up @@ -510,12 +512,13 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var
vDml_diag(i,J) = vDml(i)
enddo ; enddo

!$OMP do
!$OMP do
do j=js,je ; do k=1,nz ; do i=is,ie
h(i,j,k) = h(i,j,k) - dt*G%IareaT(i,j) * &
((uhml(I,j,k) - uhml(I-1,j,k)) + (vhml(i,J,k) - vhml(i,J-1,k)))
if (h(i,j,k) < h_min) h(i,j,k) = h_min
enddo ; enddo ; enddo
!$OMP end parallel
!$OMP end parallel

! Whenever thickness changes let the diag manager know, target grids
! for vertical remapping may need to be regenerated.
Expand Down Expand Up @@ -590,23 +593,23 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)
real :: u_star ! surface friction velocity, interpolated to velocity points [Z T-1 ~> m s-1].
real :: mom_mixrate ! rate at which momentum is homogenized within mixed layer [T-1 ~> s-1]
real :: timescale ! mixing growth timescale [T ~> s]
real :: h_min ! The minimum layer thickness [H ~> m or kg m-2]. h_min could be 0.
real :: h_neglect ! tiny thickness usually lost in roundoff and can be neglected [H ~> m or kg m-2]
real :: dz_neglect ! tiny thickness that usually lost in roundoff and can be neglected [Z ~> m]
real :: I4dt ! 1/(4 dt) [T-1 ~> s-1]
real :: I2htot ! Twice the total mixed layer thickness at velocity points [H ~> m or kg m-2]
real :: z_topx2 ! depth of the top of a layer at velocity points [H ~> m or kg m-2]
real :: hx2 ! layer thickness at velocity points [H ~> m or kg m-2]
real :: a(SZK_(GV)) ! A non-dimensional value relating the overall flux
! magnitudes (uDml & vDml) to the realized flux in a
! layer. The vertical sum of a() through the pieces of
! the mixed layer must be 0.
real :: a(SZK_(GV)) ! A non-dimensional value relating the overall flux magnitudes (uDml & vDml)
! to the realized flux in a layer [nondim]. The vertical sum of a()
! through the pieces of the mixed layer must be 0.
real :: uDml(SZIB_(G)) ! The zonal and meridional volume fluxes in the upper
real :: vDml(SZI_(G)) ! half of the mixed layer [H L2 T-1 ~> m3 s-1 or kg s-1].
real :: utimescale_diag(SZIB_(G),SZJ_(G)) ! The restratification timescales
real :: vtimescale_diag(SZI_(G),SZJB_(G)) ! in the zonal and meridional
! directions [T ~> s], stored in 2-D
real :: utimescale_diag(SZIB_(G),SZJ_(G)) ! The restratification timescales in the zonal and
real :: vtimescale_diag(SZI_(G),SZJB_(G)) ! meridional directions [T ~> s], stored in 2-D
! arrays for diagnostic purposes.
real :: uDml_diag(SZIB_(G),SZJ_(G)), vDml_diag(SZI_(G),SZJB_(G))
real :: uDml_diag(SZIB_(G),SZJ_(G)) ! A 2D copy of uDml for diagnostics [H L2 T-1 ~> m3 s-1 or kg s-1]
real :: vDml_diag(SZI_(G),SZJB_(G)) ! A 2D copy of vDml for diagnostics [H L2 T-1 ~> m3 s-1 or kg s-1]
logical :: use_EOS ! If true, density is calculated from T & S using an equation of state.

integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
Expand All @@ -619,6 +622,7 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)

if ((nkml<2) .or. (CS%ml_restrat_coef<=0.0)) return

h_min = 0.5*GV%Angstrom_H ! This should be GV%Angstrom_H, but that value would change answers.
uDml(:) = 0.0 ; vDml(:) = 0.0
I4dt = 0.25 / dt
g_Rho0 = GV%g_Earth / GV%Rho0
Expand All @@ -633,14 +637,10 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)

p0(:) = 0.0
EOSdom(:) = EOS_domain(G%HI, halo=1)
!$OMP parallel default(none) shared(is,ie,js,je,G,GV,US,htot,Rml_av,tv,p0,h,h_avail,EOSdom, &
!$OMP h_neglect,g_Rho0,I4dt,CS,uhml,uhtr,dt,vhml,vhtr, &
!$OMP utimescale_diag,vtimescale_diag,forces,dz_neglect, &
!$OMP uDml_diag,vDml_diag,nkml) &
!$OMP private(Rho0,h_vel,u_star,absf,mom_mixrate,timescale, &
!$OMP I2htot,z_topx2,hx2,a) &
!$OMP firstprivate(uDml,vDml)
!$OMP do
!$OMP parallel default(shared) private(Rho0,h_vel,u_star,absf,mom_mixrate,timescale, &
!$OMP I2htot,z_topx2,hx2,a) &
!$OMP firstprivate(uDml,vDml)
!$OMP do
do j=js-1,je+1
do i=is-1,ie+1
htot(i,j) = 0.0 ; Rml_av(i,j) = 0.0
Expand All @@ -664,7 +664,7 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)
! 2. Add exponential tail to stream-function?

! U - Component
!$OMP do
!$OMP do
do j=js,je ; do I=is-1,ie
h_vel = 0.5*(htot(i,j) + htot(i+1,j)) * GV%H_to_Z

Expand Down Expand Up @@ -711,7 +711,7 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)
enddo ; enddo

! V- component
!$OMP do
!$OMP do
do J=js-1,je ; do i=is,ie
h_vel = 0.5*(htot(i,j) + htot(i,j+1)) * GV%H_to_Z

Expand Down Expand Up @@ -756,12 +756,13 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)
vDml_diag(i,J) = vDml(i)
enddo ; enddo

!$OMP do
!$OMP do
do j=js,je ; do k=1,nkml ; do i=is,ie
h(i,j,k) = h(i,j,k) - dt*G%IareaT(i,j) * &
((uhml(I,j,k) - uhml(I-1,j,k)) + (vhml(i,J,k) - vhml(i,J-1,k)))
if (h(i,j,k) < h_min) h(i,j,k) = h_min
enddo ; enddo ; enddo
!$OMP end parallel
!$OMP end parallel

! Whenever thickness changes let the diag manager know, target grids
! for vertical remapping may need to be regenerated.
Expand Down
12 changes: 6 additions & 6 deletions src/parameterizations/vertical/MOM_diabatic_aux.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -650,7 +650,7 @@ end subroutine set_pen_shortwave


!> Diagnose a mixed layer depth (MLD) determined by a given density difference with the surface.
!> This routine is appropriate in MOM_diabatic_driver due to its position within the time stepping.
!> This routine is appropriate in MOM_diabatic_aux due to its position within the time stepping.
subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, US, diagPtr, &
id_N2subML, id_MLDsq, dz_subML)
type(ocean_grid_type), intent(in) :: G !< Grid type
Expand Down Expand Up @@ -781,7 +781,7 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, US,
end subroutine diagnoseMLDbyDensityDifference

!> Diagnose a mixed layer depth (MLD) determined by the depth a given energy value would mix.
!> This routine is appropriate in MOM_diabatic_driver due to its position within the time stepping.
!> This routine is appropriate in MOM_diabatic_aux due to its position within the time stepping.
subroutine diagnoseMLDbyEnergy(id_MLD, h, tv, G, GV, US, Mixing_Energy, diagPtr)
! Author: Brandon Reichl
! Date: October 2, 2020
Expand Down Expand Up @@ -1377,7 +1377,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
write(0,*) 'applyBoundaryFluxesInOut(): netT,netS,netH=',netHeat(i),netSalt(i),netMassInOut(i)
write(0,*) 'applyBoundaryFluxesInOut(): dT,dS,dH=',dTemp,dSalt,dThickness
write(0,*) 'applyBoundaryFluxesInOut(): h(n),h(n+1),k=',hOld,h2d(i,k),k
call MOM_error(FATAL, "MOM_diabatic_driver.F90, applyBoundaryFluxesInOut(): "//&
call MOM_error(FATAL, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//&
"Complete mass loss in column!")
endif

Expand All @@ -1392,7 +1392,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
write(0,*) 'applyBoundaryFluxesInOut(): netHeat,netSalt,netMassIn,netMassOut=',&
netHeat(i),netSalt(i),netMassIn(i),netMassOut(i)

call MOM_error(FATAL, "MOM_diabatic_driver.F90, applyBoundaryFluxesInOut(): "//&
call MOM_error(FATAL, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//&
"Mass loss over land?")
endif

Expand Down Expand Up @@ -1526,13 +1526,13 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
call forcing_SinglePointPrint(fluxes,G,iGround(i),jGround(i),'applyBoundaryFluxesInOut (grounding)')
write(mesg(1:45),'(3es15.3)') G%geoLonT( iGround(i), jGround(i) ), &
G%geoLatT( iGround(i), jGround(i)), hGrounding(i)*GV%H_to_m
call MOM_error(WARNING, "MOM_diabatic_driver.F90, applyBoundaryFluxesInOut(): "//&
call MOM_error(WARNING, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//&
"Mass created. x,y,dh= "//trim(mesg), all_print=.true.)
enddo

if (numberOfGroundings - maxGroundings > 0) then
write(mesg, '(i4)') numberOfGroundings - maxGroundings
call MOM_error(WARNING, "MOM_diabatic_driver:F90, applyBoundaryFluxesInOut(): "//&
call MOM_error(WARNING, "MOM_diabatic_aux:F90, applyBoundaryFluxesInOut(): "//&
trim(mesg) // " groundings remaining")
endif
endif
Expand Down
2 changes: 1 addition & 1 deletion src/tracer/MOM_tracer_diabatic.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -635,7 +635,7 @@ subroutine applyTracerBoundaryFluxesInOut(G, GV, Tr, dt, fluxes, h, evap_CFL_lim
if (numberOfGroundings - maxGroundings > 0) then
write(mesg, '(i4)') numberOfGroundings - maxGroundings
call MOM_error(WARNING, "MOM_tracer_vertical.F90, applyTracerBoundaryFluxesInOut(): "//&
trim(mesg) // " groundings remaining")
trim(mesg) // " groundings remaining", all_print=.true.)
endif
endif

Expand Down

0 comments on commit 4640461

Please sign in to comment.