From 5da1b8fb8e6d3efc0a6200d3786988ee0ed01e25 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 25 Jun 2021 18:02:47 -0400 Subject: [PATCH] (+)Modified some optional arguments Modified optional arguments in 4 modules to reflect their actual usage. 1. Eliminated the optional argument full_prec to zonal_flux_adjust and meridional_flux_adjust, which were always called with the hard-coded value "true", and made the optional arguments monotonic and simple_2nd to PPM_reconstruction_[xy] mandatory. 2. Eliminated the optional argument eta_bt to calculate_diagnostic_fields, which was never present. 3. Made the two optional arguments to unit_scaling_init mandatory. 4. Eliminated the optional do_i argument to F_to_ent, which was never present in calls, and made the parameter just_read_params to entrain_diffusive_init mandatory. All answers are bitwise identical. --- src/core/MOM_continuity_PPM.F90 | 157 +++++++----------- src/diagnostics/MOM_diagnostics.F90 | 13 +- src/framework/MOM_unit_scaling.F90 | 54 +++--- .../lateral/MOM_hor_visc.F90 | 3 + .../vertical/MOM_diabatic_driver.F90 | 2 - .../vertical/MOM_entrain_diffusive.F90 | 72 +++----- .../vertical/MOM_vert_friction.F90 | 2 +- 7 files changed, 117 insertions(+), 186 deletions(-) diff --git a/src/core/MOM_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90 index 7b90297c64..d8b6cddaaa 100644 --- a/src/core/MOM_continuity_PPM.F90 +++ b/src/core/MOM_continuity_PPM.F90 @@ -430,7 +430,7 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, uhbt, OBC, & if (present(uhbt)) then call zonal_flux_adjust(u, h_in, h_L, h_R, uhbt(:,j), uh_tot_0, duhdu_tot_0, du, & du_max_CFL, du_min_CFL, dt, G, GV, US, CS, visc_rem, & - j, ish, ieh, do_I, .true., uh, OBC=OBC) + j, ish, ieh, do_I, uh, OBC=OBC) if (present(u_cor)) then ; do k=1,nz do I=ish-1,ieh ; u_cor(I,j,k) = u(I,j,k) + du(I) * visc_rem(I,k) ; enddo @@ -710,7 +710,7 @@ end subroutine zonal_face_thickness !! desired barotropic (layer-summed) transport. subroutine zonal_flux_adjust(u, h_in, h_L, h_R, uhbt, uh_tot_0, duhdu_tot_0, & du, du_max_CFL, du_min_CFL, dt, G, GV, US, CS, visc_rem, & - j, ish, ieh, do_I_in, full_precision, uh_3d, OBC) + j, ish, ieh, do_I_in, uh_3d, OBC) type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure. real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1]. @@ -746,9 +746,6 @@ subroutine zonal_flux_adjust(u, h_in, h_L, h_R, uhbt, uh_tot_0, duhdu_tot_0, & integer, intent(in) :: ieh !< End of index range. logical, dimension(SZIB_(G)), intent(in) :: do_I_in !< !! A logical flag indicating which I values to work on. - logical, optional, intent(in) :: full_precision !< - !! A flag indicating how carefully to iterate. The - !! default is .true. (more accurate). real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), optional, intent(inout) :: uh_3d !< !! Volume flux through zonal faces = u*h*dy [H L2 T-1 ~> m3 s-1 or kg s-1]. type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure. @@ -768,10 +765,9 @@ subroutine zonal_flux_adjust(u, h_in, h_L, h_R, uhbt, uh_tot_0, duhdu_tot_0, & real :: tol_eta ! The tolerance for the current iteration [H ~> m or kg m-2]. real :: tol_vel ! The tolerance for velocity in the current iteration [L T-1 ~> m s-1]. integer :: i, k, nz, itt, max_itts = 20 - logical :: full_prec, domore, do_I(SZIB_(G)) + logical :: domore, do_I(SZIB_(G)) nz = GV%ke - full_prec = .true. ; if (present(full_precision)) full_prec = full_precision uh_aux(:,:) = 0.0 ; duhdu(:,:) = 0.0 @@ -787,16 +783,12 @@ subroutine zonal_flux_adjust(u, h_in, h_L, h_R, uhbt, uh_tot_0, duhdu_tot_0, & enddo do itt=1,max_itts - if (full_prec) then - select case (itt) - case (:1) ; tol_eta = 1e-6 * CS%tol_eta - case (2) ; tol_eta = 1e-4 * CS%tol_eta - case (3) ; tol_eta = 1e-2 * CS%tol_eta - case default ; tol_eta = CS%tol_eta - end select - else - tol_eta = CS%tol_eta_aux ; if (itt<=1) tol_eta = 1e-6 * CS%tol_eta_aux - endif + select case (itt) + case (:1) ; tol_eta = 1e-6 * CS%tol_eta + case (2) ; tol_eta = 1e-4 * CS%tol_eta + case (3) ; tol_eta = 1e-2 * CS%tol_eta + case default ; tol_eta = CS%tol_eta + end select tol_vel = CS%tol_vel do I=ish-1,ieh @@ -809,30 +801,23 @@ subroutine zonal_flux_adjust(u, h_in, h_L, h_R, uhbt, uh_tot_0, duhdu_tot_0, & if ((dt * min(G%IareaT(i,j),G%IareaT(i+1,j))*abs(uh_err(I)) > tol_eta) .or. & (CS%better_iter .and. ((abs(uh_err(I)) > tol_vel * duhdu_tot(I)) .or. & (abs(uh_err(I)) > uh_err_best(I))) )) then - ! Use Newton's method, provided it stays bounded. Otherwise bisect - ! the value with the appropriate bound. - if (full_prec) then - ddu = -uh_err(I) / duhdu_tot(I) - du_prev = du(I) - du(I) = du(I) + ddu - if (abs(ddu) < 1.0e-15*abs(du(I))) then - do_I(I) = .false. ! ddu is small enough to quit. - elseif (ddu > 0.0) then - if (du(I) >= du_max(I)) then - du(I) = 0.5*(du_prev + du_max(I)) - if (du_max(I) - du_prev < 1.0e-15*abs(du(I))) do_I(I) = .false. - endif - else ! ddu < 0.0 - if (du(I) <= du_min(I)) then - du(I) = 0.5*(du_prev + du_min(I)) - if (du_prev - du_min(I) < 1.0e-15*abs(du(I))) do_I(I) = .false. - endif + ! Use Newton's method, provided it stays bounded. Otherwise bisect + ! the value with the appropriate bound. + ddu = -uh_err(I) / duhdu_tot(I) + du_prev = du(I) + du(I) = du(I) + ddu + if (abs(ddu) < 1.0e-15*abs(du(I))) then + do_I(I) = .false. ! ddu is small enough to quit. + elseif (ddu > 0.0) then + if (du(I) >= du_max(I)) then + du(I) = 0.5*(du_prev + du_max(I)) + if (du_max(I) - du_prev < 1.0e-15*abs(du(I))) do_I(I) = .false. + endif + else ! ddu < 0.0 + if (du(I) <= du_min(I)) then + du(I) = 0.5*(du_prev + du_min(I)) + if (du_prev - du_min(I) < 1.0e-15*abs(du(I))) do_I(I) = .false. endif - else - ! Use Newton's method, provided it stays bounded, just like above. - du(I) = du(I) - uh_err(I) / duhdu_tot(I) - if ((du(I) >= du_max(I)) .or. (du(I) <= du_min(I))) & - du(I) = 0.5*(du_max(I) + du_min(I)) endif if (do_I(I)) domore = .true. else @@ -950,7 +935,7 @@ subroutine set_zonal_BT_cont(u, h_in, h_L, h_R, BT_cont, uh_tot_0, duhdu_tot_0, do I=ish-1,ieh ; zeros(I) = 0.0 ; enddo call zonal_flux_adjust(u, h_in, h_L, h_R, zeros, uh_tot_0, duhdu_tot_0, du0, & du_max_CFL, du_min_CFL, dt, G, GV, US, CS, visc_rem, & - j, ish, ieh, do_I, .true.) + j, ish, ieh, do_I) ! Determine the westerly- and easterly- fluxes. Choose a sufficiently ! negative velocity correction for the easterly-flux, and a sufficiently @@ -1253,7 +1238,7 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, vhbt, OBC, & if (present(vhbt)) then call meridional_flux_adjust(v, h_in, h_L, h_R, vhbt(:,J), vh_tot_0, dvhdv_tot_0, dv, & dv_max_CFL, dv_min_CFL, dt, G, GV, US, CS, visc_rem, & - j, ish, ieh, do_I, .true., vh, OBC=OBC) + j, ish, ieh, do_I, vh, OBC=OBC) if (present(v_cor)) then ; do k=1,nz do i=ish,ieh ; v_cor(i,J,k) = v(i,J,k) + dv(i) * visc_rem(i,k) ; enddo @@ -1537,7 +1522,7 @@ end subroutine merid_face_thickness !> Returns the barotropic velocity adjustment that gives the desired barotropic (layer-summed) transport. subroutine meridional_flux_adjust(v, h_in, h_L, h_R, vhbt, vh_tot_0, dvhdv_tot_0, & dv, dv_max_CFL, dv_min_CFL, dt, G, GV, US, CS, visc_rem, & - j, ish, ieh, do_I_in, full_precision, vh_3d, OBC) + j, ish, ieh, do_I_in, vh_3d, OBC) type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure. real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & @@ -1572,8 +1557,6 @@ subroutine meridional_flux_adjust(v, h_in, h_L, h_R, vhbt, vh_tot_0, dvhdv_tot_0 integer, intent(in) :: ieh !< End of index range. logical, dimension(SZI_(G)), & intent(in) :: do_I_in !< A flag indicating which I values to work on. - logical, optional, intent(in) :: full_precision !< A flag indicating how carefully to - !! iterate. The default is .true. (more accurate). real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & optional, intent(inout) :: vh_3d !< Volume flux through meridional !! faces = v*h*dx [H L2 T-1 ~> m3 s-1 or kg s-1]. @@ -1594,10 +1577,9 @@ subroutine meridional_flux_adjust(v, h_in, h_L, h_R, vhbt, vh_tot_0, dvhdv_tot_0 real :: tol_eta ! The tolerance for the current iteration [H ~> m or kg m-2]. real :: tol_vel ! The tolerance for velocity in the current iteration [L T-1 ~> m s-1]. integer :: i, k, nz, itt, max_itts = 20 - logical :: full_prec, domore, do_I(SZI_(G)) + logical :: domore, do_I(SZI_(G)) nz = GV%ke - full_prec = .true. ; if (present(full_precision)) full_prec = full_precision vh_aux(:,:) = 0.0 ; dvhdv(:,:) = 0.0 @@ -1613,16 +1595,12 @@ subroutine meridional_flux_adjust(v, h_in, h_L, h_R, vhbt, vh_tot_0, dvhdv_tot_0 enddo do itt=1,max_itts - if (full_prec) then - select case (itt) - case (:1) ; tol_eta = 1e-6 * CS%tol_eta - case (2) ; tol_eta = 1e-4 * CS%tol_eta - case (3) ; tol_eta = 1e-2 * CS%tol_eta - case default ; tol_eta = CS%tol_eta - end select - else - tol_eta = CS%tol_eta_aux ; if (itt<=1) tol_eta = 1e-6 * CS%tol_eta_aux - endif + select case (itt) + case (:1) ; tol_eta = 1e-6 * CS%tol_eta + case (2) ; tol_eta = 1e-4 * CS%tol_eta + case (3) ; tol_eta = 1e-2 * CS%tol_eta + case default ; tol_eta = CS%tol_eta + end select tol_vel = CS%tol_vel do i=ish,ieh @@ -1637,28 +1615,21 @@ subroutine meridional_flux_adjust(v, h_in, h_L, h_R, vhbt, vh_tot_0, dvhdv_tot_0 (abs(vh_err(i)) > vh_err_best(i))) )) then ! Use Newton's method, provided it stays bounded. Otherwise bisect ! the value with the appropriate bound. - if (full_prec) then - ddv = -vh_err(i) / dvhdv_tot(i) - dv_prev = dv(i) - dv(i) = dv(i) + ddv - if (abs(ddv) < 1.0e-15*abs(dv(i))) then - do_I(i) = .false. ! ddv is small enough to quit. - elseif (ddv > 0.0) then - if (dv(i) >= dv_max(i)) then - dv(i) = 0.5*(dv_prev + dv_max(i)) - if (dv_max(i) - dv_prev < 1.0e-15*abs(dv(i))) do_I(i) = .false. - endif - else ! dvv(i) < 0.0 - if (dv(i) <= dv_min(i)) then - dv(i) = 0.5*(dv_prev + dv_min(i)) - if (dv_prev - dv_min(i) < 1.0e-15*abs(dv(i))) do_I(i) = .false. - endif + ddv = -vh_err(i) / dvhdv_tot(i) + dv_prev = dv(i) + dv(i) = dv(i) + ddv + if (abs(ddv) < 1.0e-15*abs(dv(i))) then + do_I(i) = .false. ! ddv is small enough to quit. + elseif (ddv > 0.0) then + if (dv(i) >= dv_max(i)) then + dv(i) = 0.5*(dv_prev + dv_max(i)) + if (dv_max(i) - dv_prev < 1.0e-15*abs(dv(i))) do_I(i) = .false. + endif + else ! dvv(i) < 0.0 + if (dv(i) <= dv_min(i)) then + dv(i) = 0.5*(dv_prev + dv_min(i)) + if (dv_prev - dv_min(i) < 1.0e-15*abs(dv(i))) do_I(i) = .false. endif - else - ! Use Newton's method, provided it stays bounded, just like above. - dv(i) = dv(i) - vh_err(i) / dvhdv_tot(i) - if ((dv(i) >= dv_max(i)) .or. (dv(i) <= dv_min(i))) & - dv(i) = 0.5*(dv_max(i) + dv_min(i)) endif if (do_I(i)) domore = .true. else @@ -1776,7 +1747,7 @@ subroutine set_merid_BT_cont(v, h_in, h_L, h_R, BT_cont, vh_tot_0, dvhdv_tot_0, do i=ish,ieh ; zeros(i) = 0.0 ; enddo call meridional_flux_adjust(v, h_in, h_L, h_R, zeros, vh_tot_0, dvhdv_tot_0, dv0, & dv_max_CFL, dv_min_CFL, dt, G, GV, US, CS, visc_rem, & - j, ish, ieh, do_I, .true.) + j, ish, ieh, do_I) ! Determine the southerly- and northerly- fluxes. Choose a sufficiently ! negative velocity correction for the northerly-flux, and a sufficiently @@ -1871,10 +1842,10 @@ subroutine PPM_reconstruction_x(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_ type(loop_bounds_type), intent(in) :: LB !< Active loop bounds structure. real, intent(in) :: h_min !< The minimum thickness !! that can be obtained by a concave parabolic fit. - logical, optional, intent(in) :: monotonic !< If true, use the + logical, intent(in) :: monotonic !< If true, use the !! Colella & Woodward monotonic limiter. !! Otherwise use a simple positive-definite limiter. - logical, optional, intent(in) :: simple_2nd !< If true, use the + logical, intent(in) :: simple_2nd !< If true, use the !! arithmetic mean thicknesses as the default edge values !! for a simple 2nd order scheme. type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure. @@ -1884,15 +1855,11 @@ subroutine PPM_reconstruction_x(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_ real, parameter :: oneSixth = 1./6. real :: h_ip1, h_im1 real :: dMx, dMn - logical :: use_CW84, use_2nd character(len=256) :: mesg integer :: i, j, isl, iel, jsl, jel, n, stencil logical :: local_open_BC type(OBC_segment_type), pointer :: segment => NULL() - use_CW84 = .false. ; if (present(monotonic)) use_CW84 = monotonic - use_2nd = .false. ; if (present(simple_2nd)) use_2nd = simple_2nd - local_open_BC = .false. if (present(OBC)) then ; if (associated(OBC)) then local_open_BC = OBC%open_u_BCs_exist_globally @@ -1901,7 +1868,7 @@ subroutine PPM_reconstruction_x(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_ isl = LB%ish-1 ; iel = LB%ieh+1 ; jsl = LB%jsh ; jel = LB%jeh ! This is the stencil of the reconstruction, not the scheme overall. - stencil = 2 ; if (use_2nd) stencil = 1 + stencil = 2 ; if (simple_2nd) stencil = 1 if ((isl-stencil < G%isd) .or. (iel+stencil > G%ied)) then write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_x called with a ", & @@ -1916,7 +1883,7 @@ subroutine PPM_reconstruction_x(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_ call MOM_error(FATAL,mesg) endif - if (use_2nd) then + if (simple_2nd) then do j=jsl,jel ; do i=isl,iel h_im1 = G%mask2dT(i-1,j) * h_in(i-1,j) + (1.0-G%mask2dT(i-1,j)) * h_in(i,j) h_ip1 = G%mask2dT(i+1,j) * h_in(i+1,j) + (1.0-G%mask2dT(i+1,j)) * h_in(i,j) @@ -1990,7 +1957,7 @@ subroutine PPM_reconstruction_x(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_ enddo endif - if (use_CW84) then + if (monotonic) then call PPM_limit_CW84(h_in, h_L, h_R, G, isl, iel, jsl, jel) else call PPM_limit_pos(h_in, h_L, h_R, h_min, G, isl, iel, jsl, jel) @@ -2010,10 +1977,10 @@ subroutine PPM_reconstruction_y(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_ type(loop_bounds_type), intent(in) :: LB !< Active loop bounds structure. real, intent(in) :: h_min !< The minimum thickness !! that can be obtained by a concave parabolic fit. - logical, optional, intent(in) :: monotonic !< If true, use the + logical, intent(in) :: monotonic !< If true, use the !! Colella & Woodward monotonic limiter. !! Otherwise use a simple positive-definite limiter. - logical, optional, intent(in) :: simple_2nd !< If true, use the + logical, intent(in) :: simple_2nd !< If true, use the !! arithmetic mean thicknesses as the default edge values !! for a simple 2nd order scheme. type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure. @@ -2023,15 +1990,11 @@ subroutine PPM_reconstruction_y(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_ real, parameter :: oneSixth = 1./6. real :: h_jp1, h_jm1 real :: dMx, dMn - logical :: use_CW84, use_2nd character(len=256) :: mesg integer :: i, j, isl, iel, jsl, jel, n, stencil logical :: local_open_BC type(OBC_segment_type), pointer :: segment => NULL() - use_CW84 = .false. ; if (present(monotonic)) use_CW84 = monotonic - use_2nd = .false. ; if (present(simple_2nd)) use_2nd = simple_2nd - local_open_BC = .false. if (present(OBC)) then ; if (associated(OBC)) then local_open_BC = OBC%open_v_BCs_exist_globally @@ -2040,7 +2003,7 @@ subroutine PPM_reconstruction_y(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_ isl = LB%ish ; iel = LB%ieh ; jsl = LB%jsh-1 ; jel = LB%jeh+1 ! This is the stencil of the reconstruction, not the scheme overall. - stencil = 2 ; if (use_2nd) stencil = 1 + stencil = 2 ; if (simple_2nd) stencil = 1 if ((isl < G%isd) .or. (iel > G%ied)) then write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_y called with a ", & @@ -2055,7 +2018,7 @@ subroutine PPM_reconstruction_y(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_ call MOM_error(FATAL,mesg) endif - if (use_2nd) then + if (simple_2nd) then do j=jsl,jel ; do i=isl,iel h_jm1 = G%mask2dT(i,j-1) * h_in(i,j-1) + (1.0-G%mask2dT(i,j-1)) * h_in(i,j) h_jp1 = G%mask2dT(i,j+1) * h_in(i,j+1) + (1.0-G%mask2dT(i,j+1)) * h_in(i,j) @@ -2127,7 +2090,7 @@ subroutine PPM_reconstruction_y(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_ enddo endif - if (use_CW84) then + if (monotonic) then call PPM_limit_CW84(h_in, h_L, h_R, G, isl, iel, jsl, jel) else call PPM_limit_pos(h_in, h_L, h_R, h_min, G, isl, iel, jsl, jel) diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 0fbec91bc0..e6b01af33d 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -197,7 +197,7 @@ module MOM_diagnostics contains !> Diagnostics not more naturally calculated elsewhere are computed here. subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & - dt, diag_pre_sync, G, GV, US, CS, eta_bt) + dt, diag_pre_sync, G, GV, US, CS) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -227,11 +227,6 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & type(diag_grid_storage), intent(in) :: diag_pre_sync !< Target grids from previous timestep type(diagnostics_CS), intent(inout) :: CS !< Control structure returned by a !! previous call to diagnostics_init. - real, dimension(SZI_(G),SZJ_(G)), & - optional, intent(in) :: eta_bt !< An optional barotropic - !! variable that gives the "correct" free surface height (Boussinesq) or total water column - !! mass per unit area (non-Boussinesq). This is used to dilate the layer thicknesses when - !! calculating interface heights [H ~> m or kg m-2]. ! Local variables real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: usq ! squared eastward velocity [L2 T-2 ~> m2 s-2] @@ -390,7 +385,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & endif if (associated(CS%e)) then - call find_eta(h, tv, G, GV, US, CS%e, eta_bt) + call find_eta(h, tv, G, GV, US, CS%e) if (CS%id_e > 0) call post_data(CS%id_e, CS%e, CS%diag) endif @@ -400,7 +395,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & CS%e_D(i,j,k) = CS%e(i,j,k) + G%bathyT(i,j) enddo ; enddo ; enddo else - call find_eta(h, tv, G, GV, US, CS%e_D, eta_bt) + call find_eta(h, tv, G, GV, US, CS%e_D) do k=1,nz+1 ; do j=js,je ; do i=is,ie CS%e_D(i,j,k) = CS%e_D(i,j,k) + G%bathyT(i,j) enddo ; enddo ; enddo @@ -1935,7 +1930,7 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag call wave_speed_init(CS%wave_speed_CSp, remap_answers_2018=remap_answers_2018, & better_speed_est=better_speed_est, min_speed=wave_speed_min, & wave_speed_tol=wave_speed_tol) - call wave_speed_init(CS%wave_speed_CSp, remap_answers_2018=remap_answers_2018) +!### call wave_speed_init(CS%wave_speed_CSp, remap_answers_2018=remap_answers_2018) call safe_alloc_ptr(CS%cg1,isd,ied,jsd,jed) if (CS%id_Rd1>0) call safe_alloc_ptr(CS%Rd1,isd,ied,jsd,jed) if (CS%id_Rd_ebt>0) call safe_alloc_ptr(CS%Rd1,isd,ied,jsd,jed) diff --git a/src/framework/MOM_unit_scaling.F90 b/src/framework/MOM_unit_scaling.F90 index fea1ac4910..dbcd2405ec 100644 --- a/src/framework/MOM_unit_scaling.F90 +++ b/src/framework/MOM_unit_scaling.F90 @@ -54,8 +54,8 @@ module MOM_unit_scaling !> Allocates and initializes the ocean model unit scaling type subroutine unit_scaling_init( param_file, US ) - type(param_file_type), optional, intent(in) :: param_file !< Parameter file handle/type - type(unit_scale_type), optional, pointer :: US !< A dimensional unit scaling type + type(param_file_type), intent(in) :: param_file !< Parameter file handle/type + type(unit_scale_type), pointer :: US !< A dimensional unit scaling type ! This routine initializes a unit_scale_type structure (US). @@ -66,39 +66,33 @@ subroutine unit_scaling_init( param_file, US ) # include "version_variable.h" character(len=16) :: mdl = "MOM_unit_scaling" - if (.not.present(US)) return - if (associated(US)) call MOM_error(FATAL, & 'unit_scaling_init: called with an associated US pointer.') allocate(US) - if (present(param_file)) then - ! Read all relevant parameters and write them to the model log. - call log_version(param_file, mdl, version, & - "Parameters for doing unit scaling of variables.", debugging=.true.) - call get_param(param_file, mdl, "Z_RESCALE_POWER", Z_power, & - "An integer power of 2 that is used to rescale the model's "//& - "internal units of depths and heights. Valid values range from -300 to 300.", & - units="nondim", default=0, debuggingParam=.true.) - call get_param(param_file, mdl, "L_RESCALE_POWER", L_power, & - "An integer power of 2 that is used to rescale the model's "//& - "internal units of lateral distances. Valid values range from -300 to 300.", & - units="nondim", default=0, debuggingParam=.true.) - call get_param(param_file, mdl, "T_RESCALE_POWER", T_power, & - "An integer power of 2 that is used to rescale the model's "//& - "internal units of time. Valid values range from -300 to 300.", & - units="nondim", default=0, debuggingParam=.true.) - call get_param(param_file, mdl, "R_RESCALE_POWER", R_power, & - "An integer power of 2 that is used to rescale the model's "//& - "internal units of density. Valid values range from -300 to 300.", & + ! Read all relevant parameters and write them to the model log. + call log_version(param_file, mdl, version, & + "Parameters for doing unit scaling of variables.", debugging=.true.) + call get_param(param_file, mdl, "Z_RESCALE_POWER", Z_power, & + "An integer power of 2 that is used to rescale the model's "//& + "internal units of depths and heights. Valid values range from -300 to 300.", & + units="nondim", default=0, debuggingParam=.true.) + call get_param(param_file, mdl, "L_RESCALE_POWER", L_power, & + "An integer power of 2 that is used to rescale the model's "//& + "internal units of lateral distances. Valid values range from -300 to 300.", & + units="nondim", default=0, debuggingParam=.true.) + call get_param(param_file, mdl, "T_RESCALE_POWER", T_power, & + "An integer power of 2 that is used to rescale the model's "//& + "internal units of time. Valid values range from -300 to 300.", & + units="nondim", default=0, debuggingParam=.true.) + call get_param(param_file, mdl, "R_RESCALE_POWER", R_power, & + "An integer power of 2 that is used to rescale the model's "//& + "internal units of density. Valid values range from -300 to 300.", & + units="nondim", default=0, debuggingParam=.true.) + call get_param(param_file, mdl, "Q_RESCALE_POWER", Q_power, & + "An integer power of 2 that is used to rescale the model's "//& + "internal units of heat content. Valid values range from -300 to 300.", & units="nondim", default=0, debuggingParam=.true.) - call get_param(param_file, mdl, "Q_RESCALE_POWER", Q_power, & - "An integer power of 2 that is used to rescale the model's "//& - "internal units of heat content. Valid values range from -300 to 300.", & - units="nondim", default=0, debuggingParam=.true.) - else - Z_power = 0 ; L_power = 0 ; T_power = 0 ; R_power = 0 ; Q_power = 0 - endif if (abs(Z_power) > 300) call MOM_error(FATAL, "unit_scaling_init: "//& "Z_RESCALE_POWER is outside of the valid range of -300 to 300.") diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 24afcf8cd8..c588a1faa4 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -2509,6 +2509,7 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, MEKE, ADp) if (CS%Laplacian .or. get_all) then endif end subroutine hor_visc_init + !> Calculates factors in the anisotropic orientation tensor to be align with the grid. !! With n1=1 and n2=0, this recovers the approach of Large et al, 2001. subroutine align_aniso_tensor_to_grid(CS, n1, n2) @@ -2525,6 +2526,7 @@ subroutine align_aniso_tensor_to_grid(CS, n1, n2) CS%n1n1_m_n2n2_h(:,:) = ( n1 * n1 - n2 * n2 ) * recip_n2_norm CS%n1n1_m_n2n2_q(:,:) = ( n1 * n1 - n2 * n2 ) * recip_n2_norm end subroutine align_aniso_tensor_to_grid + !> Apply a 1-1-4-1-1 Laplacian filter one time on GME diffusive flux to reduce any !! horizontal two-grid-point noise subroutine smooth_GME(CS,G,GME_flux_h,GME_flux_q) @@ -2589,6 +2591,7 @@ subroutine smooth_GME(CS,G,GME_flux_h,GME_flux_q) endif enddo ! s-loop end subroutine smooth_GME + !> Deallocates any variables allocated in hor_visc_init. subroutine hor_visc_end(CS) type(hor_visc_CS), pointer :: CS !< The control structure returned by a diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index e579661f47..7a75802a84 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -3084,8 +3084,6 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di 'Salinity', 'PSU') endif - - !call set_diffusivity_init(Time, G, param_file, diag, CS%set_diff_CSp, CS%int_tide_CSp) CS%id_Kd_int = register_diag_field('ocean_model', 'Kd_interface', diag%axesTi, Time, & 'Total diapycnal diffusivity at interfaces', 'm2 s-1', conversion=US%Z2_T_to_m2_s) if (CS%use_energetic_PBL) then diff --git a/src/parameterizations/vertical/MOM_entrain_diffusive.F90 b/src/parameterizations/vertical/MOM_entrain_diffusive.F90 index a558f9dd2b..32cdce4d2a 100644 --- a/src/parameterizations/vertical/MOM_entrain_diffusive.F90 +++ b/src/parameterizations/vertical/MOM_entrain_diffusive.F90 @@ -352,7 +352,7 @@ subroutine entrainment_diffusive(h, tv, fluxes, dt, G, GV, US, CS, ea, eb, & ! what maxF(kb+1) should be. do i=is,ie ; min_eakb(i) = MIN(htot(i), max_eakb(i)) ; enddo call find_maxF_kb(h_bl, Sref, Ent_bl, I_dSkbp1, min_eakb, max_eakb, & - kmb, is, ie, G, GV, CS, F_kb_maxEnt, do_i_in = do_i) + kmb, is, ie, G, GV, CS, F_kb_maxEnt, do_i_in=do_i) do i=is,ie do_entrain_eakb = .false. @@ -891,7 +891,7 @@ end subroutine entrainment_diffusive !> This subroutine calculates the actual entrainments (ea and eb) and the !! amount of surface forcing that is applied to each layer if there is no bulk !! mixed layer. -subroutine F_to_ent(F, h, kb, kmb, j, G, GV, CS, dsp1_ds, eakb, Ent_bl, ea, eb, do_i_in) +subroutine F_to_ent(F, h, kb, kmb, j, G, GV, CS, dsp1_ds, eakb, Ent_bl, ea, eb) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure real, dimension(SZI_(G),SZK_(GV)), intent(in) :: F !< The density flux through a layer within @@ -920,31 +920,13 @@ subroutine F_to_ent(F, h, kb, kmb, j, G, GV, CS, dsp1_ds, eakb, Ent_bl, ea, eb, real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(inout) :: eb !< The amount of fluid entrained from the layer !! below within this time step [H ~> m or kg m-2]. - logical, dimension(SZI_(G)), & - optional, intent(in) :: do_i_in !< Indicates which i-points to work on. -! This subroutine calculates the actual entrainments (ea and eb) and the -! amount of surface forcing that is applied to each layer if there is no bulk -! mixed layer. real :: h1 ! The thickness in excess of the minimum that will remain ! after exchange with the layer below [H ~> m or kg m-2]. - logical :: do_i(SZI_(G)) integer :: i, k, is, ie, nz is = G%isc ; ie = G%iec ; nz = GV%ke - if (present(do_i_in)) then - do i=is,ie ; do_i(i) = do_i_in(i) ; enddo - do i=G%isc,G%iec ; if (do_i(i)) then - is = i ; exit - endif ; enddo - do i=G%iec,G%isc,-1 ; if (do_i(i)) then - ie = i ; exit - endif ; enddo - else - do i=is,ie ; do_i(i) = .true. ; enddo - endif - do i=is,ie ea(i,j,nz) = 0.0 ; eb(i,j,nz) = 0.0 enddo @@ -952,7 +934,7 @@ subroutine F_to_ent(F, h, kb, kmb, j, G, GV, CS, dsp1_ds, eakb, Ent_bl, ea, eb, do i=is,ie eb(i,j,kmb) = max(2.0*Ent_bl(i,Kmb+1) - eakb(i), 0.0) enddo - do k=nz-1,kmb+1,-1 ; do i=is,ie ; if (do_i(i)) then + do k=nz-1,kmb+1,-1 ; do i=is,ie if (k > kb(i)) then ! With a bulk mixed layer, surface buoyancy fluxes are applied ! elsewhere, so F should always be nonnegative. @@ -970,9 +952,9 @@ subroutine F_to_ent(F, h, kb, kmb, j, G, GV, CS, dsp1_ds, eakb, Ent_bl, ea, eb, ! up into the buffer layer. eb(i,j,k) = eb(i,j,k+1) + max(0.0, h(i,j,k+1) - GV%Angstrom_H) endif - endif ; enddo ; enddo + enddo ; enddo k = kmb - do i=is,ie ; if (do_i(i)) then + do i=is,ie ! Adjust the previously calculated entrainment from below by the deepest ! buffer layer to account for entrainment of thin interior layers . if (kb(i) > kmb+1) & @@ -981,8 +963,8 @@ subroutine F_to_ent(F, h, kb, kmb, j, G, GV, CS, dsp1_ds, eakb, Ent_bl, ea, eb, ! Determine the entrainment from above for each buffer layer. h1 = (h(i,j,k) - GV%Angstrom_H) + (eb(i,j,k) - ea(i,j,k+1)) ea(i,j,k) = MAX(Ent_bl(i,K), Ent_bl(i,K)-0.5*h1, -h1) - endif ; enddo - do k=kmb-1,2,-1 ; do i=is,ie ; if (do_i(i)) then + enddo + do k=kmb-1,2,-1 ; do i=is,ie ! Determine the entrainment from below for each buffer layer. eb(i,j,k) = max(2.0*Ent_bl(i,K+1) - ea(i,j,k+1), 0.0) @@ -992,11 +974,11 @@ subroutine F_to_ent(F, h, kb, kmb, j, G, GV, CS, dsp1_ds, eakb, Ent_bl, ea, eb, ! if (h1 >= 0.0) then ; ea(i,j,k) = Ent_bl(i,K) ! elseif (Ent_bl(i,K)+0.5*h1 >= 0.0) then ; ea(i,j,k) = Ent_bl(i,K)-0.5*h1 ! else ; ea(i,j,k) = -h1 ; endif - endif ; enddo ; enddo - do i=is,ie ; if (do_i(i)) then + enddo ; enddo + do i=is,ie eb(i,j,1) = max(2.0*Ent_bl(i,2) - ea(i,j,2), 0.0) ea(i,j,1) = 0.0 - endif ; enddo + enddo else ! not BULKMIXEDLAYER ! Calculate the entrainment by each layer from above and below. ! Entrainment is always positive, but F may be negative due to @@ -1511,7 +1493,7 @@ subroutine F_kb_to_ea_kb(h_bl, Sref, Ent_bl, I_dSkbp1, F_kb, kmb, i, & ! the maximum. zeros(i) = 0.0 call find_maxF_kb(h_bl, Sref, Ent_bl, I_dSkbp1, zeros, ea_kb, & - kmb, i, i, G, GV, CS, maxF, ent_maxF, F_thresh = F_kb) + kmb, i, i, G, GV, CS, maxF, ent_maxF, F_thresh=F_kb) err_max = dS_kbp1 * maxF(i) - val ! If err_max is negative, there is no good solution, so use the maximum ! value of F in the valid range. @@ -1693,7 +1675,7 @@ subroutine determine_Ea_kb(h_bl, dtKd_kb, Sref, I_dSkbp1, Ent_bl, ea_kbp1, & do_any = .false. ; do i=is,ie ; if (redo_i(i)) do_any = .true. ; enddo if (.not.do_any) exit call determine_dSkb(h_bl, Sref, Ent_bl, Ent, is, ie, kmb, G, GV, .true., dS_kb, & - ddSkb_dE, dS_lay, ddSlay_dE, do_i_in = redo_i) + ddSkb_dE, dS_lay, ddSlay_dE, do_i_in=redo_i) do i=is,ie ; if (redo_i(i)) then ! The correct root is bracketed between E_min and E_max. ! Note the following limits: Ent >= 0 ; fa > 1 ; fk > 0 @@ -1757,7 +1739,7 @@ subroutine determine_Ea_kb(h_bl, dtKd_kb, Sref, I_dSkbp1, Ent_bl, ea_kbp1, & ! Update the value of dS_kb for consistency with Ent. if (present(F_kb) .or. present(dFdfm_kb)) & call determine_dSkb(h_bl, Sref, Ent_bl, Ent, is, ie, kmb, G, GV, .true., & - dS_kb, do_i_in = do_i) + dS_kb, do_i_in=do_i) if (present(F_kb)) then ; do i=is,ie ; if (do_i(i)) then F_kb(i) = Ent(i) * (dS_kb(i) * I_dSkbp1(i)) @@ -1878,7 +1860,7 @@ subroutine find_maxF_kb(h_bl, Sref, Ent_bl, I_dSkbp1, min_ent_in, max_ent_in, & do i=ie1,is,-1 ; if (do_i(i)) is1 = i ; enddo ! Find the value of F and its derivative at min_ent. call determine_dSkb(h_bl, Sref, Ent_bl, minent, is1, ie1, kmb, G, GV, .false., & - dS_kb, ddSkb_dE, do_i_in = do_i) + dS_kb, ddSkb_dE, do_i_in=do_i) do i=is1,ie1 ; if (do_i(i)) then F_minent(i) = minent(i) * dS_kb(i) * I_dSkbp1(i) dF_dE_min(i) = (dS_kb(i) + minent(i)*ddSkb_dE(i)) * I_dSkbp1(i) @@ -1958,7 +1940,7 @@ subroutine find_maxF_kb(h_bl, Sref, Ent_bl, I_dSkbp1, min_ent_in, max_ent_in, & endif call determine_dSkb(h_bl, Sref, Ent_bl, ent, is1, ie1, kmb, G, GV, .false., & - dS_kb, ddSkb_dE, do_i_in = do_i) + dS_kb, ddSkb_dE, do_i_in=do_i) do i=is1,ie1 ; if (do_i(i)) then F(i) = ent(i)*dS_kb(i)*I_dSkbp1(i) dF_dent(i) = (dS_kb(i) + ent(i)*ddSkb_dE(i)) * I_dSkbp1(i) @@ -2050,8 +2032,7 @@ subroutine find_maxF_kb(h_bl, Sref, Ent_bl, I_dSkbp1, min_ent_in, max_ent_in, & enddo if (doany) then ! For efficiency, could save previous value of dS_anom_lim_best? - call determine_dSkb(h_bl, Sref, Ent_bl, ent_best, is, ie, kmb, G, GV, .true., & - dS_kb_lim) + call determine_dSkb(h_bl, Sref, Ent_bl, ent_best, is, ie, kmb, G, GV, .true., dS_kb_lim) do i=is,ie F_best(i) = ent_best(i)*dS_kb_lim(i)*I_dSkbp1(i) ! The second test seems necessary because of roundoff differences that @@ -2088,14 +2069,13 @@ subroutine entrain_diffusive_init(Time, G, GV, US, param_file, diag, CS, just_re !! output. type(entrain_diffusive_CS), pointer :: CS !< A pointer that is set to point to the control !! structure. - logical, optional, intent(in) :: just_read_params !< If present and true, this call will - !! only read parameters logging them or registering + logical, intent(in) :: just_read_params !< If true, this call will only read + !! and log parameters without registering !! any diagnostics ! Local variables real :: dt ! The dynamics timestep, used here in the default for TOLERANCE_ENT, in MKS units [s] real :: Kd ! A diffusivity used in the default for TOLERANCE_ENT, in MKS units [m2 s-1] - logical :: just_read ! If true, just read parameters but do nothing else. ! This include declares and sets the variable "version". # include "version_variable.h" character(len=40) :: mdl = "MOM_entrain_diffusive" ! This module's name. @@ -2107,30 +2087,28 @@ subroutine entrain_diffusive_init(Time, G, GV, US, param_file, diag, CS, just_re endif allocate(CS) - just_read = .false. ; if (present(just_read_params)) just_read = just_read_params - CS%diag => diag CS%bulkmixedlayer = (GV%nkml > 0) -! Set default, read and log parameters - if (.not.just_read) call log_version(param_file, mdl, version, "") + ! Set default, read and log parameters + if (.not.just_read_params) call log_version(param_file, mdl, version, "") call get_param(param_file, mdl, "MAX_ENT_IT", CS%max_ent_it, & "The maximum number of iterations that may be used to "//& - "calculate the interior diapycnal entrainment.", default=5, do_not_log=just_read) + "calculate the interior diapycnal entrainment.", default=5, do_not_log=just_read_params) ! In this module, KD is only used to set the default for TOLERANCE_ENT. [m2 s-1] call get_param(param_file, mdl, "KD", Kd, default=0.0) call get_param(param_file, mdl, "DT", dt, & "The (baroclinic) dynamics time step.", units = "s", & - fail_if_missing=.true., do_not_log=just_read) + fail_if_missing=.true., do_not_log=just_read_params) call get_param(param_file, mdl, "TOLERANCE_ENT", CS%Tolerance_Ent, & "The tolerance with which to solve for entrainment values.", & units="m", default=MAX(100.0*GV%Angstrom_m,1.0e-4*sqrt(dt*Kd)), scale=GV%m_to_H, & - do_not_log=just_read) + do_not_log=just_read_params) CS%Rho_sig_off = 1000.0*US%kg_m3_to_R - if (.not.just_read) then + if (.not.just_read_params) then CS%id_Kd = register_diag_field('ocean_model', 'Kd_effective', diag%axesTL, Time, & 'Diapycnal diffusivity as applied', 'm2 s-1', conversion=US%Z2_T_to_m2_s) CS%id_diff_work = register_diag_field('ocean_model', 'diff_work', diag%axesTi, Time, & @@ -2138,7 +2116,7 @@ subroutine entrain_diffusive_init(Time, G, GV, US, param_file, diag, CS, just_re 'W m-2', conversion=US%RZ3_T3_to_W_m2) endif - if (just_read) deallocate(CS) + if (just_read_params) deallocate(CS) end subroutine entrain_diffusive_init diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 6465df1d4e..5b85c5f5f6 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -1888,7 +1888,7 @@ end subroutine vertvisc_init subroutine updateCFLtruncationValue(Time, CS, activate) type(time_type), target, intent(in) :: Time !< Current model time type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure - logical, optional, intent(in) :: activate !< Specifiy whether to record the value of + logical, optional, intent(in) :: activate !< Specify whether to record the value of !! Time as the beginning of the ramp period ! Local variables