Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

(*)Avoid negative thicknesses in mixed_layer_restrat #76

Merged
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
85 changes: 43 additions & 42 deletions src/parameterizations/lateral/MOM_mixed_layer_restrat.F90
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
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
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