From 5da1b8fb8e6d3efc0a6200d3786988ee0ed01e25 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 25 Jun 2021 18:02:47 -0400 Subject: [PATCH 1/8] (+)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 From 52fe57631f4aa33251a024fe8718fa0959befc02 Mon Sep 17 00:00:00 2001 From: He Wang <35150900+herrwang0@users.noreply.github.com> Date: Tue, 29 Jun 2021 14:59:45 -0400 Subject: [PATCH 2/8] Bugfix: Minor changes on clipping topography from file (#1428) * Minor changes to topography initialized from files This is intended to extend the compatibility with negative topography. 1. A bug that occurs when clipping topography with a negative MINIMUM_DEPTH is fixed. 2. MASKING_DEPTH can now be negative. 3. A warning will be given if MASKING_DEPTH is set to be smaller than MINIMUM_DEPTH. * Adding an exception to avoid answer changes To keep answer unchnaged in the test cases (at land points), D is clipped at 0.5*min_depth when min_depth > 0. * Change the if-statement on whether MASK_DEPTH is default to a simple comparison * Change the if-statement on whether MASK_DEPTH is default * Change the if-statement comparison on MASK_DEPTH Reverted back to a simple comparison. A mistake in a previous commit in MOM_grid_initialize.F90 is corrected. Co-authored-by: Robert Hallberg --- src/initialization/MOM_grid_initialize.F90 | 12 ++++++--- .../MOM_shared_initialization.F90 | 26 +++++++++++++------ 2 files changed, 27 insertions(+), 11 deletions(-) diff --git a/src/initialization/MOM_grid_initialize.F90 b/src/initialization/MOM_grid_initialize.F90 index 0fac3e15b4..fbee77d130 100644 --- a/src/initialization/MOM_grid_initialize.F90 +++ b/src/initialization/MOM_grid_initialize.F90 @@ -9,7 +9,7 @@ module MOM_grid_initialize use MOM_domains, only : To_North, To_South, To_East, To_West use MOM_domains, only : MOM_domain_type, clone_MOM_domain, deallocate_MOM_domain use MOM_dyn_horgrid, only : dyn_horgrid_type, set_derived_dyn_horgrid -use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, is_root_pe +use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave use MOM_file_parser, only : get_param, log_param, log_version, param_file_type use MOM_io, only : MOM_read_data, slasher, file_exists, stdout @@ -1217,11 +1217,17 @@ subroutine initialize_masks(G, PF, US) units="m", default=0.0, scale=m_to_Z_scale) call get_param(PF, mdl, "MASKING_DEPTH", mask_depth, & "The depth below which to mask points as land points, for which all "//& - "fluxes are zeroed out. MASKING_DEPTH is ignored if negative.", & + "fluxes are zeroed out. MASKING_DEPTH needs to be smaller than MINIMUM_DEPTH", & units="m", default=-9999.0, scale=m_to_Z_scale) + if (mask_depth > min_depth) then + mask_depth = -9999.0*m_to_Z_scale + call MOM_error(WARNING, "MOM_grid_init: initialize_masks "//& + 'MASKING_DEPTH is larger than MINIMUM_DEPTH and therefore ignored.') + endif + Dmin = min_depth - if (mask_depth>=0.) Dmin = mask_depth + if (mask_depth /= -9999.*m_to_Z_scale) Dmin = mask_depth G%mask2dCu(:,:) = 0.0 ; G%mask2dCv(:,:) = 0.0 ; G%mask2dBu(:,:) = 0.0 diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90 index daaefd4b98..336a85d5bc 100644 --- a/src/initialization/MOM_shared_initialization.F90 +++ b/src/initialization/MOM_shared_initialization.F90 @@ -413,17 +413,27 @@ subroutine limit_topography(D, G, param_file, max_depth, US) "The depth below which to mask the ocean as land.", & units="m", default=-9999.0, scale=m_to_Z, do_not_log=.true.) -! Make sure that min_depth < D(x,y) < max_depth - if (mask_depth < -9990.*m_to_Z) then - do j=G%jsd,G%jed ; do i=G%isd,G%ied - D(i,j) = min( max( D(i,j), 0.5*min_depth ), max_depth ) - enddo ; enddo + if (mask_depth > min_depth) then + mask_depth = -9999.0*m_to_Z + call MOM_error(WARNING, "MOM_shared_initialization: limit_topography "//& + 'MASKING_DEPTH is larger than MINIMUM_DEPTH and therefore ignored.') + endif + + ! Make sure that min_depth < D(x,y) < max_depth for ocean points + if (mask_depth == -9999.*m_to_Z) then + if (min_depth > 0.0) then ! This is retained to avoid answer changes (over the land points) in the test cases. + do j=G%jsd,G%jed ; do i=G%isd,G%ied + D(i,j) = min( max( D(i,j), 0.5*min_depth ), max_depth ) + enddo ; enddo + else + do j=G%jsd,G%jed ; do i=G%isd,G%ied + D(i,j) = min( max( D(i,j), min_depth ), max_depth ) + enddo ; enddo + endif else do j=G%jsd,G%jed ; do i=G%isd,G%ied - if (D(i,j)>0.) then + if (D(i,j) > mask_depth) then D(i,j) = min( max( D(i,j), min_depth ), max_depth ) - else - D(i,j) = 0. endif enddo ; enddo endif From 4e3cc5d286bb5c6cfde9b6995b5c16610c67ff92 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Tue, 29 Jun 2021 23:00:08 -0400 Subject: [PATCH 3/8] Update gitlab pipeline to use scripts from MOM6-examples - Reduce entries in .gitlab-ci.yml to mostly just one line by invoking scripts in MOM6-examples - This undoes early direction of moving detailed control into the pipeline yaml - todo: define a library build process from within each child repository --- .gitlab-ci.yml | 183 +++++++++++++++++-------------------------------- 1 file changed, 64 insertions(+), 119 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 5d84c0c176..5e1da1c1f9 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -1,113 +1,78 @@ stages: - - merge+setup - builds - run - tests - cleanup +variables: + CACHE_DIR: "/lustre/f2/scratch/oar.gfdl.ogrp-account/runner/cache/" + + # Merges MOM6 with dev/gfdl. Changes directory to test directory, if it exists. +# - set cache location +# - get MOM6-examples/tools/MRS scripts by cloning Gaea-stats and then MOM6-examples +# - set working directory to MOM6-examples +# - pull down latest of dev/gfdl (MOM6-examples might be ahead of Gaea-stats) before_script: - - MOM6_SRC=$CI_PROJECT_DIR - - echo Cache directory set to ${CACHE_DIR:=/lustre/f2/scratch/oar.gfdl.ogrp-account/runner/cache/} - - git pull --no-edit https://github.com/NOAA-GFDL/MOM6.git dev/gfdl && git submodule init && git submodule update - - pwd ; ls + - echo Cache directory set to $CACHE_DIR + - echo -e "\e[0Ksection_start:`date +%s`:before[collapsed=true]\r\e[0KPre-script" + - git clone https://gitlab.gfdl.noaa.gov/ogrp/Gaea-stats-MOM6-examples.git tests + - cd tests && git submodule init && git submodule update + - cd MOM6-examples && git checkout dev/gfdl && git pull + - echo -e "\e[0Ksection_end:`date +%s`:before\r\e[0K" # Tests that merge with dev/gfdl works. merge: - stage: merge+setup + stage: builds tags: - ncrc4 script: - - pwd ; ls + - cd $CI_PROJECT_DIR - git pull --no-edit https://github.com/NOAA-GFDL/MOM6.git dev/gfdl -# Clones regression repo, if necessary, pulls latest of everything, and sets up working space -setup: - stage: merge+setup - tags: - - ncrc4 - script: - - pwd ; ls - # Clone regressions directory - - git clone --recursive http://gitlab.gfdl.noaa.gov/ogrp/Gaea-stats-MOM6-examples.git tests && cd tests - # Install / update testing scripts - - git clone -b new-code-struct https://github.com/adcroft/MRS.git MRS - # Update MOM6-examples and submodules - - (cd MOM6-examples && git checkout . && git checkout dev/gfdl && git pull && git submodule init && git submodule update) - - (cd MOM6-examples/src/MOM6 && git submodule update) - - test -d MOM6-examples/src/LM3 || make -f MRS/Makefile.clone clone_gfdl -s - - make -f MRS/Makefile.clone MOM6-examples/.datasets -s - - env > gitlab_session.log - # Show hashes for final setup - - git show --oneline - - git submodule status - - (cd MOM6-examples && git submodule status --recursive src) - # Cache everything under tests to unpack for each subsequent stage - - cd ../ ; time tar zcf $CACHE_DIR/tests_$CI_PIPELINE_ID.tgz tests - # Compiles gnu:repro: stage: builds tags: - ncrc4 script: - - time tar zxf $CACHE_DIR/tests_$CI_PIPELINE_ID.tgz && cd tests - - time make -f MRS/Makefile.build MOM6_SRC=../ build_gnu -s -j - - time make -f MRS/Makefile.build MOM6_SRC=../ static_gnu -s -j - - time tar zvcf $CACHE_DIR/build-gnu-repro-$CI_PIPELINE_ID.tgz `find build/gnu -name MOM6` + - time make -f tools/MRS/Makefile MOM6_SRC=../.. pipeline-build-repro-gnu -s -j + - time make -f tools/MRS/Makefile MOM6_SRC=../.. pipeline-build-static-gnu -s -j gnu:ocean-only-nolibs: stage: builds tags: - ncrc4 script: - - time tar zxf $CACHE_DIR/tests_$CI_PIPELINE_ID.tgz && cd tests - - make -f MRS/Makefile.build build/gnu/env && cd build/gnu - # mkdir -p build/gnu/repro/symmetric_dynamic/ocean_only && cd build/gnu/repro/symmetric_dynamic/ocean_only - - ../../MOM6-examples/src/mkmf/bin/list_paths -l ../../../config_src/{drivers/solo_driver,memory/dynamic_symmetric,infra/FMS1,ext*} ../../../src ../../MOM6-examples/src/FMS - - sed -i '/FMS\/.*\/test_/d' path_names - - ../../MOM6-examples/src/mkmf/bin/mkmf -t ../../MOM6-examples/src/mkmf/templates/ncrc-gnu.mk -p MOM6 -c"-Duse_libMPI -Duse_netCDF" path_names - - time (source ./env ; make NETCDF=3 REPRO=1 MOM6 -s -j) + - make -f tools/MRS/Makefile pipeline-build-gnu-oceanonly-nolibs gnu:ice-ocean-nolibs: stage: builds tags: - ncrc4 script: - - time tar zxf $CACHE_DIR/tests_$CI_PIPELINE_ID.tgz && cd tests - - make -f MRS/Makefile.build build/gnu/env && cd build/gnu - # mkdir -p build/gnu/repro/symmetric_dynamic/ocean_only && cd build/gnu/repro/symmetric_dynamic/ocean_only - - ../../MOM6-examples/src/mkmf/bin/list_paths -l ../../../config_src/{drivers/FMS_cap,memory/dynamic_nonsymmetric,infra/FMS1,ext*} ../../../src ../../MOM6-examples/src/{FMS,coupler,SIS2,icebergs,ice_param,land_null,atmos_null} - - sed -i '/FMS\/.*\/test_/d' path_names - - ../../MOM6-examples/src/mkmf/bin/mkmf -t ../../MOM6-examples/src/mkmf/templates/ncrc-gnu.mk -p MOM6 -c"-Duse_libMPI -Duse_netCDF -D_USE_LEGACY_LAND_ -Duse_AM3_physics" path_names - - time (source ./env ; make NETCDF=3 REPRO=1 MOM6 -s -j) + - make -f tools/MRS/Makefile pipeline-build-gnu-iceocean-nolibs intel:repro: stage: builds tags: - ncrc4 script: - - time tar zxf $CACHE_DIR/tests_$CI_PIPELINE_ID.tgz && cd tests - - make -f MRS/Makefile.build MOM6_SRC=../ build_intel -s -j - - time tar zvcf $CACHE_DIR/build-intel-repro-$CI_PIPELINE_ID.tgz `find build/intel -name MOM6` + - time make -f tools/MRS/Makefile MOM6_SRC=../.. pipeline-build-repro-intel -s -j pgi:repro: stage: builds tags: - ncrc4 script: - - time tar zxf $CACHE_DIR/tests_$CI_PIPELINE_ID.tgz && cd tests - - make -f MRS/Makefile.build MOM6_SRC=../ build_pgi -s -j - - time tar zvcf $CACHE_DIR/build-pgi-repro-$CI_PIPELINE_ID.tgz `find build/pgi -name MOM6` + - time make -f tools/MRS/Makefile MOM6_SRC=../.. pipeline-build-repro-pgi -s -j gnu:debug: stage: builds tags: - ncrc4 script: - - time tar zxf $CACHE_DIR/tests_$CI_PIPELINE_ID.tgz && cd tests - - make -f MRS/Makefile.build MOM6_SRC=../ debug_gnu -s -j - - time tar zvcf $CACHE_DIR/build-gnu-debug-$CI_PIPELINE_ID.tgz `find build/gnu -name MOM6` + - time make -f tools/MRS/Makefile MOM6_SRC=../.. pipeline-build-debug-gnu -s -j # Runs run: @@ -115,41 +80,43 @@ run: tags: - ncrc4 script: - - time tar zxf $CACHE_DIR/tests_$CI_PIPELINE_ID.tgz && cd tests - - time tar zxf $CACHE_DIR/build-gnu-repro-$CI_PIPELINE_ID.tgz - - time tar zxf $CACHE_DIR/build-intel-repro-$CI_PIPELINE_ID.tgz - - time tar zxf $CACHE_DIR/build-pgi-repro-$CI_PIPELINE_ID.tgz - # time tar zxf $CACHE_DIR/build-gnu-debug-$CI_PIPELINE_ID.tgz - - (echo '#!/bin/tcsh';echo 'make -f MRS/Makefile.tests all') > job.sh - - sbatch --clusters=c3,c4 --nodes=29 --time=0:34:00 --account=gfdl_o --qos=debug --job-name=mom6_regressions --output=log.$CI_PIPELINE_ID --wait job.sh || MJOB_RETURN_STATE=Fail - - cat log.$CI_PIPELINE_ID - - test -z "$MJOB_RETURN_STATE" - - test -f restart_results_gnu.tar.gz - - time tar zvcf $CACHE_DIR/results-$CI_PIPELINE_ID.tgz *.tar.gz + - make -f tools/MRS/Makefile mom6-pipeline-run gnu.testing: stage: run tags: - ncrc4 + before_script: + - echo -e "\e[0Ksection_start:`date +%s`:submodules[collapsed=true]\r\e[0KCloning submodules" + - git submodule init ; git submodule update + - echo -e "\e[0Ksection_end:`date +%s`:submodules\r\e[0K" script: + - echo -e "\e[0Ksection_start:`date +%s`:compile[collapsed=true]\r\e[0KCompiling executables" - cd .testing - module unload PrgEnv-pgi PrgEnv-intel PrgEnv-gnu darshan ; module load PrgEnv-gnu ; module unload netcdf gcc ; module load gcc/7.3.0 cray-hdf5 cray-netcdf - make work/local-env - make -s -j + - echo -e "\e[0Ksection_end:`date +%s`:compile\r\e[0K" - (echo '#!/bin/bash';echo '. ./work/local-env/bin/activate';echo 'make MPIRUN="srun -mblock --exclusive" test -s -j') > job.sh - - sbatch --clusters=c3,c4 --nodes=5 --time=0:05:00 --account=gfdl_o --qos=debug --job-name=MOM6.gnu.testing --output=log.$CI_PIPELINE_ID --wait job.sh || cat log.$CI_PIPELINE_ID && make test + - sbatch --clusters=c3,c4 --nodes=5 --time=0:05:00 --account=gfdl_o --qos=debug --job-name=MOM6.gnu.testing --output=log.$CI_PIPELINE_ID --wait job.sh && make test || cat log.$CI_PIPELINE_ID intel.testing: stage: run tags: - ncrc4 + before_script: + - echo -e "\e[0Ksection_start:`date +%s`:submodules[collapsed=true]\r\e[0KCloning submodules" + - git submodule init ; git submodule update + - echo -e "\e[0Ksection_end:`date +%s`:submodules\r\e[0K" script: + - echo -e "\e[0Ksection_start:`date +%s`:compile[collapsed=true]\r\e[0KCompiling executables" - cd .testing - module unload PrgEnv-pgi PrgEnv-intel PrgEnv-gnu darshan; module load PrgEnv-intel; module unload netcdf intel; module load intel/18.0.6.288 cray-hdf5 cray-netcdf - make work/local-env - make -s -j + - echo -e "\e[0Ksection_end:`date +%s`:compile\r\e[0K" - (echo '#!/bin/bash';echo '. ./work/local-env/bin/activate';echo 'make MPIRUN="srun -mblock --exclusive" test -s -j') > job.sh - - sbatch --clusters=c3,c4 --nodes=5 --time=0:05:00 --account=gfdl_o --qos=debug --job-name=MOM6.gnu.testing --output=log.$CI_PIPELINE_ID --wait job.sh || cat log.$CI_PIPELINE_ID && make test + - sbatch --clusters=c3,c4 --nodes=5 --time=0:05:00 --account=gfdl_o --qos=debug --job-name=MOM6.gnu.testing --output=log.$CI_PIPELINE_ID --wait job.sh && make test || cat log.$CI_PIPELINE_ID # Tests gnu:non-symmetric: @@ -157,113 +124,91 @@ gnu:non-symmetric: tags: - ncrc4 script: - - time tar zxf $CACHE_DIR/tests_$CI_PIPELINE_ID.tgz && cd tests - - time tar zxf $CACHE_DIR/results-$CI_PIPELINE_ID.tgz - - make -f MRS/Makefile.tests gnu_non_symmetric + - make -f tools/MRS/Makefile mom6-pipeline-test-gnu_non_symmetric -intel:non-symmetric: +gnu:symmetric: stage: tests tags: - ncrc4 script: - - time tar zxf $CACHE_DIR/tests_$CI_PIPELINE_ID.tgz && cd tests - - time tar zxf $CACHE_DIR/results-$CI_PIPELINE_ID.tgz - - make -f MRS/Makefile.tests intel_non_symmetric + - make -f tools/MRS/Makefile mom6-pipeline-test-gnu_symmetric -pgi:non-symmetric: +gnu:memory: stage: tests tags: - ncrc4 script: - - time tar zxf $CACHE_DIR/tests_$CI_PIPELINE_ID.tgz && cd tests - - time tar zxf $CACHE_DIR/results-$CI_PIPELINE_ID.tgz - - make -f MRS/Makefile.tests pgi_non_symmetric + - make -f tools/MRS/Makefile mom6-pipeline-test-gnu_memory -gnu:symmetric: +gnu:static: stage: tests tags: - ncrc4 script: - - time tar zxf $CACHE_DIR/tests_$CI_PIPELINE_ID.tgz && cd tests - - time tar zxf $CACHE_DIR/results-$CI_PIPELINE_ID.tgz - - make -f MRS/Makefile.tests gnu_symmetric + - make -f tools/MRS/Makefile mom6-pipeline-test-gnu_static -intel:symmetric: +gnu:restart: stage: tests tags: - ncrc4 script: - - time tar zxf $CACHE_DIR/tests_$CI_PIPELINE_ID.tgz && cd tests - - time tar zxf $CACHE_DIR/results-$CI_PIPELINE_ID.tgz - - make -f MRS/Makefile.tests intel_symmetric + - make -f tools/MRS/Makefile mom6-pipeline-test-gnu_restarts -pgi:symmetric: +gnu:params: stage: tests tags: - ncrc4 script: - - time tar zxf $CACHE_DIR/tests_$CI_PIPELINE_ID.tgz && cd tests - - time tar zxf $CACHE_DIR/results-$CI_PIPELINE_ID.tgz - - make -f MRS/Makefile.tests pgi_symmetric + - make -f tools/MRS/Makefile mom6-pipeline-test-params_gnu_symmetric + allow_failure: true -gnu:layout: +intel:symmetric: stage: tests tags: - ncrc4 script: - - time tar zxf $CACHE_DIR/tests_$CI_PIPELINE_ID.tgz && cd tests - - time tar zxf $CACHE_DIR/results-$CI_PIPELINE_ID.tgz - - make -f MRS/Makefile.tests gnu_layout + - make -f tools/MRS/Makefile mom6-pipeline-test-intel_symmetric -intel:layout: +intel:non-symmetric: stage: tests tags: - ncrc4 script: - - time tar zxf $CACHE_DIR/tests_$CI_PIPELINE_ID.tgz && cd tests - - time tar zxf $CACHE_DIR/results-$CI_PIPELINE_ID.tgz - - make -f MRS/Makefile.tests intel_layout + - make -f tools/MRS/Makefile mom6-pipeline-test-intel_non_symmetric -pgi:layout: +intel:memory: stage: tests tags: - ncrc4 script: - - time tar zxf $CACHE_DIR/tests_$CI_PIPELINE_ID.tgz && cd tests - - time tar zxf $CACHE_DIR/results-$CI_PIPELINE_ID.tgz - - make -f MRS/Makefile.tests pgi_layout + - make -f tools/MRS/Makefile mom6-pipeline-test-intel_memory -gnu:static: +pgi:symmetric: stage: tests tags: - ncrc4 script: - - time tar zxf $CACHE_DIR/tests_$CI_PIPELINE_ID.tgz && cd tests - - time tar zxf $CACHE_DIR/results-$CI_PIPELINE_ID.tgz - - make -f MRS/Makefile.tests gnu_static + - make -f tools/MRS/Makefile mom6-pipeline-test-pgi_symmetric -gnu:restart: +pgi:non-symmetric: stage: tests tags: - ncrc4 script: - - time tar zxf $CACHE_DIR/tests_$CI_PIPELINE_ID.tgz && cd tests - - time tar zxf $CACHE_DIR/results-$CI_PIPELINE_ID.tgz - - make -f MRS/Makefile.tests gnu_check_restarts + - make -f tools/MRS/Makefile mom6-pipeline-test-pgi_non_symmetric -gnu:params: +pgi:memory: stage: tests tags: - ncrc4 script: - - time tar zxf $CACHE_DIR/tests_$CI_PIPELINE_ID.tgz && cd tests - - time tar zxf $CACHE_DIR/results-$CI_PIPELINE_ID.tgz - - make -f MRS/Makefile.tests params_gnu_symmetric - allow_failure: true + - make -f tools/MRS/Makefile mom6-pipeline-test-pgi_memory cleanup: stage: cleanup tags: - ncrc4 + before_script: + - echo Skipping submodule update script: - rm $CACHE_DIR/*$CI_PIPELINE_ID.tgz From d9bdbc3431abdc2fe36d75071fd9ef0ec43c2a06 Mon Sep 17 00:00:00 2001 From: Keith Lindsay Date: Fri, 9 Jul 2021 12:45:06 -0600 Subject: [PATCH 4/8] fix dim_names assignment in MOM_io:read_var_sizes fixes #1439 --- src/framework/MOM_io.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90 index 3ad2c92f41..f8cfb09382 100644 --- a/src/framework/MOM_io.F90 +++ b/src/framework/MOM_io.F90 @@ -714,7 +714,7 @@ subroutine read_var_sizes(filename, varname, ndims, sizes, match_case, caller, d if (status /= NF90_NOERR) call MOM_error(WARNING, trim(hdr) // trim(NF90_STRERROR(status)) //& " Getting dimension length for "//trim(varname)//" in "//trim(filename)) if (present(dim_names)) then - if (n <= size(dim_names)) dim_names = trim(dimname) + if (n <= size(dim_names)) dim_names(n) = trim(dimname) endif enddo deallocate(dimids) From 7c27bfba2bf27c6d78bd13792906c4438731c186 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 9 Jul 2021 18:09:11 -0400 Subject: [PATCH 5/8] +Add query_wave_properties & fix NUOPC wave queries Added a new routine, query_wave_properties, that can be call to get information about the wave properties from the waves control structure, and added a get_param call for SURFBAND_WAVENUMBERS to MOM_wave_interface_init when it is using the options that are typical with NUOPC coupler. These changes should allow the NUOPC coupler to compile and work again, while still keeping the same level of opacity in the wave_parameters_CS. All answers should be bitwise identical, although the order of some entries in the MOM_parameter_doc files may change when coupled with waves is WAVE_METHOD=SURFACE_BANDS and SURFBAND_SOURCE=COUPLER. --- config_src/drivers/nuopc_cap/mom_cap.F90 | 8 ++++-- .../nuopc_cap/mom_ocean_model_nuopc.F90 | 6 ---- src/user/MOM_wave_interface.F90 | 28 +++++++++++++++++++ 3 files changed, 33 insertions(+), 9 deletions(-) diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index 2d79674606..6561f63fd1 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -28,8 +28,9 @@ module MOM_cap_mod use MOM_get_input, only: get_MOM_input, directories use MOM_domains, only: pass_var use MOM_error_handler, only: MOM_error, FATAL, is_root_pe -use MOM_ocean_model_nuopc, only: ice_ocean_boundary_type use MOM_grid, only: ocean_grid_type, get_global_grid_size +use MOM_wave_interface, only: query_wave_properties +use MOM_ocean_model_nuopc, only: ice_ocean_boundary_type use MOM_ocean_model_nuopc, only: ocean_model_restart, ocean_public_type, ocean_state_type use MOM_ocean_model_nuopc, only: ocean_model_init_sfc use MOM_ocean_model_nuopc, only: ocean_model_init, update_ocean_model, ocean_model_end @@ -696,7 +697,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) Ice_ocean_boundary%frunoff = 0.0 if (ocean_state%use_waves) then - Ice_ocean_boundary%num_stk_bands=ocean_state%Waves%NumBands + call query_wave_properties(ocean_state%Waves, NumBands=Ice_ocean_boundary%num_stk_bands) allocate ( Ice_ocean_boundary% ustk0 (isc:iec,jsc:jec), & Ice_ocean_boundary% vstk0 (isc:iec,jsc:jec), & Ice_ocean_boundary% ustkb (isc:iec,jsc:jec,Ice_ocean_boundary%num_stk_bands), & @@ -704,7 +705,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) Ice_ocean_boundary%stk_wavenumbers (Ice_ocean_boundary%num_stk_bands)) Ice_ocean_boundary%ustk0 = 0.0 Ice_ocean_boundary%vstk0 = 0.0 - Ice_ocean_boundary%stk_wavenumbers = ocean_state%Waves%WaveNum_Cen + call query_wave_properties(ocean_state%Waves, WaveNumbers=Ice_ocean_boundary%stk_wavenumbers, & + US=ocean_state%US) Ice_ocean_boundary%ustkb = 0.0 Ice_ocean_boundary%vstkb = 0.0 endif diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 index 1d5de0dd3e..06fa905a23 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 @@ -391,12 +391,6 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i ! MOM_wave_interface_init is called regardless of the value of USE_WAVES because ! it also initializes statistical waves. call MOM_wave_interface_init(OS%Time, OS%grid, OS%GV, OS%US, param_file, OS%Waves, OS%diag) - if (OS%use_waves) then - ! I do not know why this is being set here. It seems out of place. -RWH - call get_param(param_file,mdl,"SURFBAND_WAVENUMBERS", OS%Waves%WaveNum_Cen, & - "Central wavenumbers for surface Stokes drift bands.", & - units='rad/m', default=0.12566, scale=OS%US%Z_to_m) - endif if (associated(OS%grid%Domain%maskmap)) then call initialize_ocean_public_type(OS%grid%Domain%mpp_domain, Ocean_sfc, & diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90 index e1e4ab0f77..d51f6ae46a 100644 --- a/src/user/MOM_wave_interface.F90 +++ b/src/user/MOM_wave_interface.F90 @@ -24,6 +24,7 @@ module MOM_wave_interface #include public MOM_wave_interface_init ! Public interface to fully initialize the wave routines. +public query_wave_properties ! Public interface to obtain information from the waves control structure. public Update_Surface_Waves ! Public interface to update wave information at the ! coupler/driver level. public Update_Stokes_Drift ! Public interface to update the Stokes drift profiles @@ -322,6 +323,9 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag ) "STK_BAND_COUPLER is the number of Stokes drift bands in the coupler. "// & "This has to be consistent with the number of Stokes drift bands in WW3, "//& "or the model will fail.",units='', default=1) + call get_param(param_file, mdl, "SURFBAND_WAVENUMBERS", CS%WaveNum_Cen, & + "Central wavenumbers for surface Stokes drift bands.", & + units='rad/m', default=0.12566, scale=US%Z_to_m) allocate( CS%WaveNum_Cen(CS%NumBands) ) allocate( CS%STKx0(G%isdB:G%iedB,G%jsd:G%jed,CS%NumBands)) allocate( CS%STKy0(G%isdB:G%iedB,G%jsd:G%jed,CS%NumBands)) @@ -428,6 +432,30 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag ) end subroutine MOM_wave_interface_init +!> This interface provides the caller with information from the waves control structure. +subroutine query_wave_properties(CS, NumBands, WaveNumbers, US) + type(wave_parameters_CS), pointer :: CS !< Wave parameter Control structure + integer, optional, intent(out) :: NumBands !< If present, this returns the number of + !!< wavenumber partitions in the wave discretization + real, dimension(:), optional, intent(out) :: Wavenumbers !< If present this returns the characteristic + !! wavenumbers of the wave discretization [m-1 or Z-1 ~> m-1] + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type that is used to undo + !! the dimensional scaling of the output variables, if present + integer :: n + + if (present(NumBands)) NumBands = CS%NumBands + if (present(Wavenumbers)) then + if (size(Wavenumbers) < CS%NumBands) call MOM_error(FATAL, "query_wave_properties called "//& + "with a Wavenumbers array that is smaller than the number of bands.") + if (present(US)) then + do n=1,CS%NumBands ; Wavenumbers(n) = US%m_to_Z * CS%WaveNum_Cen(n) ; enddo + else + do n=1,CS%NumBands ; Wavenumbers(n) = CS%WaveNum_Cen(n) ; enddo + endif + endif + +end subroutine query_wave_properties + !> Subroutine that handles updating of surface wave/Stokes drift related properties subroutine Update_Surface_Waves(G, GV, US, Day, dt, CS, forces) type(wave_parameters_CS), pointer :: CS !< Wave parameter Control structure From f6c1fc73a700d26b4fc84014ad69b93f956524ec Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 12 Jul 2021 18:47:14 +0000 Subject: [PATCH 6/8] Allocate Waves%WaveNum_Cen before reading into it Moved the recently added call to read SURFBAND_WAVENUMBERS into an array when SURFBAND_SOURCE=COUPLER down by a line to follow the allocate call for that array. This change should avert a memory access problem that would otherwise arise when exercising this newly added code. All answers are bitwise identical in any case that ran with the previous version, and they should reproduce answers from before this PR as a whole, although obviously this code is not as well tested as would be ideal, based on the fact that this bug was in the last commit, which also passed testing. --- src/user/MOM_wave_interface.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90 index d51f6ae46a..b612ef4270 100644 --- a/src/user/MOM_wave_interface.F90 +++ b/src/user/MOM_wave_interface.F90 @@ -323,10 +323,10 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag ) "STK_BAND_COUPLER is the number of Stokes drift bands in the coupler. "// & "This has to be consistent with the number of Stokes drift bands in WW3, "//& "or the model will fail.",units='', default=1) + allocate( CS%WaveNum_Cen(CS%NumBands) ) call get_param(param_file, mdl, "SURFBAND_WAVENUMBERS", CS%WaveNum_Cen, & "Central wavenumbers for surface Stokes drift bands.", & units='rad/m', default=0.12566, scale=US%Z_to_m) - allocate( CS%WaveNum_Cen(CS%NumBands) ) allocate( CS%STKx0(G%isdB:G%iedB,G%jsd:G%jed,CS%NumBands)) allocate( CS%STKy0(G%isdB:G%iedB,G%jsd:G%jed,CS%NumBands)) CS%WaveNum_Cen(:) = 0.0 From 01b51e825a8312b0730394f565641f4e38189953 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 15 Jul 2021 10:23:37 -0400 Subject: [PATCH 7/8] +Add query_ocean_state Added the new routine query_ocean_state to mom_ocean_model_nuopc to allow for the wave properties to be obtained by the mom_cap without having to rely on the elements of an otherwise opaque type being public. I believe that all elements of the ocean_state_type could now be declared as private, and that the version of MOM6 with the NUOPC coupler should now work as intended. --- config_src/drivers/nuopc_cap/mom_cap.F90 | 16 +++++----- .../nuopc_cap/mom_ocean_model_nuopc.F90 | 29 +++++++++++++++++-- 2 files changed, 36 insertions(+), 9 deletions(-) diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index 6561f63fd1..394cf05285 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -29,12 +29,11 @@ module MOM_cap_mod use MOM_domains, only: pass_var use MOM_error_handler, only: MOM_error, FATAL, is_root_pe use MOM_grid, only: ocean_grid_type, get_global_grid_size -use MOM_wave_interface, only: query_wave_properties use MOM_ocean_model_nuopc, only: ice_ocean_boundary_type use MOM_ocean_model_nuopc, only: ocean_model_restart, ocean_public_type, ocean_state_type use MOM_ocean_model_nuopc, only: ocean_model_init_sfc use MOM_ocean_model_nuopc, only: ocean_model_init, update_ocean_model, ocean_model_end -use MOM_ocean_model_nuopc, only: get_ocean_grid, get_eps_omesh +use MOM_ocean_model_nuopc, only: get_ocean_grid, get_eps_omesh, query_ocean_state use MOM_cap_time, only: AlarmInit use MOM_cap_methods, only: mom_import, mom_export, mom_set_geomtype, mod2med_areacor use MOM_cap_methods, only: med2mod_areacor, state_diagnose @@ -422,6 +421,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) character(len=64) :: logmsg logical :: isPresent, isPresentDiro, isPresentLogfile, isSet logical :: existflag + logical :: use_waves ! If true, the wave modules are active. integer :: userRc integer :: localPet integer :: localPeCount @@ -696,8 +696,9 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) Ice_ocean_boundary%lrunoff = 0.0 Ice_ocean_boundary%frunoff = 0.0 - if (ocean_state%use_waves) then - call query_wave_properties(ocean_state%Waves, NumBands=Ice_ocean_boundary%num_stk_bands) + call query_ocean_state(ocean_state, use_waves=use_waves) + if (use_waves) then + call query_ocean_state(ocean_state, NumWaveBands=Ice_ocean_boundary%num_stk_bands) allocate ( Ice_ocean_boundary% ustk0 (isc:iec,jsc:jec), & Ice_ocean_boundary% vstk0 (isc:iec,jsc:jec), & Ice_ocean_boundary% ustkb (isc:iec,jsc:jec,Ice_ocean_boundary%num_stk_bands), & @@ -705,11 +706,12 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) Ice_ocean_boundary%stk_wavenumbers (Ice_ocean_boundary%num_stk_bands)) Ice_ocean_boundary%ustk0 = 0.0 Ice_ocean_boundary%vstk0 = 0.0 - call query_wave_properties(ocean_state%Waves, WaveNumbers=Ice_ocean_boundary%stk_wavenumbers, & - US=ocean_state%US) + call query_ocean_state(ocean_state, WaveNumbers=Ice_ocean_boundary%stk_wavenumbers, unscale=.true.) Ice_ocean_boundary%ustkb = 0.0 Ice_ocean_boundary%vstkb = 0.0 endif + ! Consider adding this: + ! if (.not.use_waves) Ice_ocean_boundary%num_stk_bands = 0 ocean_internalstate%ptr%ocean_state_type_ptr => ocean_state call ESMF_GridCompSetInternalState(gcomp, ocean_internalstate, rc) @@ -754,7 +756,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) !These are not currently used and changing requires a nuopc dictionary change !call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_runoff_heat_flx" , "will provide") !call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_calving_heat_flx" , "will provide") - if (ocean_state%use_waves) then + if (use_waves) then if (Ice_ocean_boundary%num_stk_bands > 3) then call MOM_error(FATAL, "Number of Stokes Bands > 3, NUOPC cap not set up for this") endif diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 index 06fa905a23..03f2b1ed9d 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 @@ -58,7 +58,7 @@ module MOM_ocean_model_nuopc use mpp_mod, only : mpp_chksum use MOM_EOS, only : gsw_sp_from_sr, gsw_pt_from_ct use MOM_wave_interface, only : wave_parameters_CS, MOM_wave_interface_init -use MOM_wave_interface, only : Update_Surface_Waves +use MOM_wave_interface, only : Update_Surface_Waves, query_wave_properties use MOM_surface_forcing_nuopc, only : surface_forcing_init, convert_IOB_to_fluxes use MOM_surface_forcing_nuopc, only : convert_IOB_to_forces, ice_ocn_bnd_type_chksum use MOM_surface_forcing_nuopc, only : ice_ocean_boundary_type, surface_forcing_CS @@ -80,7 +80,7 @@ module MOM_ocean_model_nuopc public ocean_model_restart public ice_ocn_bnd_type_chksum public ocean_public_type_chksum -public get_ocean_grid +public get_ocean_grid, query_ocean_state public get_eps_omesh !> This type is used for communication with other components via the FMS coupler. @@ -1000,6 +1000,31 @@ subroutine ocean_model_flux_init(OS, verbosity) end subroutine ocean_model_flux_init +!> This interface allows certain properties that are stored in the ocean_state_type to be +!! obtained. +subroutine query_ocean_state(OS, use_waves, NumWaveBands, Wavenumbers, unscale) + type(ocean_state_type), intent(in) :: OS !< The structure with the complete ocean state + logical, optional, intent(out) :: use_waves !< Indicates whether surface waves are in use + integer, optional, intent(out) :: NumWaveBands !< If present, this gives the number of + !! wavenumber partitions in the wave discretization + real, dimension(:), optional, intent(out) :: Wavenumbers !< If present, this gives the characteristic + !! wavenumbers of the wave discretization [m-1 or Z-1 ~> m-1] + logical, optional, intent(in) :: unscale !< If present and true, undo any dimensional + !! rescaling and return dimensional values in MKS units + + logical :: undo_scaling + undo_scaling = .false. ; if (present(unscale)) undo_scaling = unscale + + if (present(use_waves)) use_waves = OS%use_waves + if (present(NumWaveBands)) call query_wave_properties(OS%Waves, NumBands=NumWaveBands) + if (present(Wavenumbers) .and. undo_scaling) then + call query_wave_properties(OS%Waves, WaveNumbers=WaveNumbers, US=OS%US) + elseif (present(Wavenumbers)) then + call query_wave_properties(OS%Waves, WaveNumbers=WaveNumbers) + endif + +end subroutine query_ocean_state + !> Ocean_stock_pe - returns the integrated stocks of heat, water, etc. for conservation checks. !! Because of the way FMS is coded, only the root PE has the integrated amount, !! while all other PEs get 0. From 2ff2419eb086a6eac8d02322ea2dabd734de556f Mon Sep 17 00:00:00 2001 From: jiandewang Date: Sun, 18 Jul 2021 13:03:38 -0400 Subject: [PATCH 8/8] initialize CS%WaveNum_Cen before read in this parameter --- src/user/MOM_wave_interface.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90 index b612ef4270..38aa6b13a5 100644 --- a/src/user/MOM_wave_interface.F90 +++ b/src/user/MOM_wave_interface.F90 @@ -324,15 +324,15 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag ) "This has to be consistent with the number of Stokes drift bands in WW3, "//& "or the model will fail.",units='', default=1) allocate( CS%WaveNum_Cen(CS%NumBands) ) - call get_param(param_file, mdl, "SURFBAND_WAVENUMBERS", CS%WaveNum_Cen, & - "Central wavenumbers for surface Stokes drift bands.", & - units='rad/m', default=0.12566, scale=US%Z_to_m) allocate( CS%STKx0(G%isdB:G%iedB,G%jsd:G%jed,CS%NumBands)) allocate( CS%STKy0(G%isdB:G%iedB,G%jsd:G%jed,CS%NumBands)) CS%WaveNum_Cen(:) = 0.0 CS%STKx0(:,:,:) = 0.0 CS%STKy0(:,:,:) = 0.0 CS%PartitionMode = 0 + call get_param(param_file, mdl, "SURFBAND_WAVENUMBERS", CS%WaveNum_Cen, & + "Central wavenumbers for surface Stokes drift bands.", & + units='rad/m', default=0.12566, scale=US%Z_to_m) case (INPUT_STRING)! A method to input the Stokes band (globally uniform) CS%DataSource = INPUT call get_param(param_file,mdl,"SURFBAND_NB",CS%NumBands, &