From 464046138500854a95a0e7f89ac1677e71c8462a Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 24 Feb 2022 13:51:55 -0500 Subject: [PATCH] (*)Avoid negative thicknesses in mixed_layer_restrat 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. --- .../lateral/MOM_mixed_layer_restrat.F90 | 85 ++++++++++--------- .../vertical/MOM_diabatic_aux.F90 | 12 +-- src/tracer/MOM_tracer_diabatic.F90 | 2 +- 3 files changed, 50 insertions(+), 49 deletions(-) diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index 04982d7171..6176f4aa1d 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -156,15 +156,16 @@ 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 @@ -172,21 +173,25 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var 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) ) @@ -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.) & @@ -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 @@ -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))) @@ -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))) @@ -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. @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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. diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 13d25f06f5..6f1988c295 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 diff --git a/src/tracer/MOM_tracer_diabatic.F90 b/src/tracer/MOM_tracer_diabatic.F90 index c865e645ad..66adfb3695 100644 --- a/src/tracer/MOM_tracer_diabatic.F90 +++ b/src/tracer/MOM_tracer_diabatic.F90 @@ -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