From a3d62f021659fdbe510cf41d871136884c3bd492 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Fri, 20 Nov 2020 18:24:38 -0800 Subject: [PATCH 01/12] Autoconf: netCDF flag fixes Autoconf netCDF flag configuration was broken due to using nc-config for Fortran flags, and assuming that C and Fortran files (modules, libraries, etc) were installed in the same directory. This patch replaces the lower level --includedir and --libdir flags with the higher level --fflags and --flibs flags, which provide more processed Fortran output. One significant issue here is that --flibs combines the directory (-L) and library (-l) flags into a single output, even though these should be split when passed to `ld`, with autoconf passing the former to LDFLAGS and the latter to LIBS. We resolve this by stripping the -l flags (via sed) and pass the remaining flags (almost always -L) to LDFLAGS. We rely on autoconf macros to assign the -l flags to LIBS. This should resolve issues where nc-fortran and nf-fortran are installed in separate directories. This was detected in a spack-based environment, where this would be commonplace. --- ac/configure.ac | 20 ++++++++++++++------ ac/deps/configure.fms.ac | 20 ++++++++++++++------ 2 files changed, 28 insertions(+), 12 deletions(-) diff --git a/ac/configure.ac b/ac/configure.ac index 14249fc09a..b069fdd56f 100644 --- a/ac/configure.ac +++ b/ac/configure.ac @@ -82,12 +82,20 @@ AX_FC_CHECK_MODULE([mpi], # netCDF configuration -AC_PATH_PROG([NC_CONFIG], [nc-config]) -AS_IF([test -n "$NC_CONFIG"], - [CPPFLAGS="$CPPFLAGS -I$($NC_CONFIG --includedir)" - FCFLAGS="$FCFLAGS -I$($NC_CONFIG --includedir)" - LDFLAGS="$LDFLAGS -L$($NC_CONFIG --libdir)"], - [AC_MSG_ERROR([Could not find nc-config.])]) + +# NOTE: `nf-config --flibs` combines library paths (-L) and libraries (-l), +# even though these ought to be separated in the invocation of `ld`. +# +# We use `sed` to strip the -l and pass the -L to LDFLAGS, and rely on autoconf +# to configure the -l flags. +AC_PROG_SED + +AC_PATH_PROG([NF_CONFIG], [nf-config]) +AS_IF([test -n "$NF_CONFIG"], + [CPPFLAGS="$CPPFLAGS $($NF_CONFIG --fflags)" + FCFLAGS="$FCFLAGS $($NF_CONFIG --fflags)" + LDFLAGS="$LDFLAGS $($NF_CONFIG --flibs | $SED -e 's/-l[[^ ]]*//g')"], + [AC_MSG_ERROR([Could not find nf-config.])]) AX_FC_CHECK_MODULE([netcdf], [], [AC_MSG_ERROR([Could not find netcdf module.])]) diff --git a/ac/deps/configure.fms.ac b/ac/deps/configure.fms.ac index 714f20df6b..4ac86f6445 100644 --- a/ac/deps/configure.fms.ac +++ b/ac/deps/configure.fms.ac @@ -73,12 +73,20 @@ AC_DEFINE([use_libMPI]) # netCDF configuration -AC_PATH_PROG([NC_CONFIG], [nc-config]) -AS_IF([test -n "$NC_CONFIG"], - [CPPFLAGS="$CPPFLAGS -I$($NC_CONFIG --includedir)" - FCFLAGS="$FCFLAGS -I$($NC_CONFIG --includedir)" - LDFLAGS="$LDFLAGS -L$($NC_CONFIG --libdir)"], - [AC_MSG_ERROR([Could not find nc-config.])]) + +# NOTE: `nf-config --flibs` combines library paths (-L) and libraries (-l), +# even though these ought to be separated in the invocation of `ld`. +# +# We use `sed` to strip the -l and pass the -L to LDFLAGS, and rely on autoconf +# to configure the -l flags. +AC_PROG_SED + +AC_PATH_PROG([NF_CONFIG], [nf-config]) +AS_IF([test -n "$NF_CONFIG"], + [CPPFLAGS="$CPPFLAGS $($NF_CONFIG --fflags)" + FCFLAGS="$FCFLAGS $($NF_CONFIG --fflags)" + LDFLAGS="$LDFLAGS $($NF_CONFIG --flibs | $SED -e 's/-l[[^ ]]*//g')"], + [AC_MSG_ERROR([Could not find nf-config.])]) AX_FC_CHECK_MODULE([netcdf], [], [AC_MSG_ERROR([Could not find netcdf module.])]) From 75983903fa690eb1c18dde8dbd58e988b21d4cb3 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 24 Nov 2020 19:36:41 -0500 Subject: [PATCH 02/12] +Add triDiagTS_Eulerian & tracer_vertdiff_Eulerian Added the new subroutines triDiagTS_Eulerian and tracer_vertdiff_Eulerian to handle the cases where there is tracer mixing but no net transfer of mass with one less 3-d array argument. Also added the new optional argument zero_mix to find_uv_at_h to reproduce the previous answers with no mixing in some cases where arrays of zeros were passed in as arguments. The calls in diabatic_ALE and diabatic_ALE_legacy were modified to use these new interfaces. All answers are bitwise identical, but there are new publicly visible routines and a new optional arguments to a publicly visible routine. --- .../vertical/MOM_diabatic_aux.F90 | 78 +++++- .../vertical/MOM_diabatic_driver.F90 | 39 ++- src/tracer/MOM_tracer_diabatic.F90 | 227 ++++++++++++++++-- 3 files changed, 300 insertions(+), 44 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index fbc6bb2d16..40e7f16e76 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -30,7 +30,7 @@ module MOM_diabatic_aux #include public diabatic_aux_init, diabatic_aux_end -public make_frazil, adjust_salt, differential_diffuse_T_S, triDiagTS +public make_frazil, adjust_salt, differential_diffuse_T_S, triDiagTS, triDiagTS_Eulerian public find_uv_at_h, diagnoseMLDbyDensityDifference, applyBoundaryFluxesInOut, set_pen_shortwave public diagnoseMLDbyEnergy @@ -395,9 +395,10 @@ subroutine triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, T, S) real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: S !< Layer salinities [ppt]. ! Local variables - real :: b1(SZIB_(G)), d1(SZIB_(G)) ! b1, c1, and d1 are variables used by the - real :: c1(SZIB_(G),SZK_(G)) ! tridiagonal solver. - real :: h_tr, b_denom_1 + real :: b1(SZIB_(G)) ! A variable used by the tridiagonal solver [H-1 ~> m-2 or m2 kg-1]. + real :: d1(SZIB_(G)) ! A variable used by the tridiagonal solver [nondim]. + real :: c1(SZIB_(G),SZK_(G)) ! A variable used by the tridiagonal solver [nondim]. + real :: h_tr, b_denom_1 ! Two temporary thicknesses [H ~> m or kg m-2]. integer :: i, j, k !$OMP parallel do default(shared) private(h_tr,b1,d1,c1,b_denom_1) @@ -425,9 +426,58 @@ subroutine triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, T, S) enddo end subroutine triDiagTS +!> This is a simple tri-diagonal solver for T and S, with mixing across interfaces but no net +!! transfer of mass. +subroutine triDiagTS_Eulerian(G, GV, is, ie, js, je, hold, ent, T, S) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + integer, intent(in) :: is !< The start i-index to work on. + integer, intent(in) :: ie !< The end i-index to work on. + integer, intent(in) :: js !< The start j-index to work on. + integer, intent(in) :: je !< The end j-index to work on. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: hold !< The layer thicknesses before entrainment, + !! [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in) :: ent !< The amount of fluid mixed across an interface + !! within this time step [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: T !< Layer potential temperatures [degC]. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: S !< Layer salinities [ppt]. + + ! Local variables + real :: b1(SZIB_(G)) ! A variable used by the tridiagonal solver [H-1 ~> m-2 or m2 kg-1]. + real :: d1(SZIB_(G)) ! A variable used by the tridiagonal solver [nondim]. + real :: c1(SZIB_(G),SZK_(G)) ! A variable used by the tridiagonal solver [nondim]. + real :: h_tr, b_denom_1 ! Two temporary thicknesses [H ~> m or kg m-2]. + integer :: i, j, k + + !$OMP parallel do default(shared) private(h_tr,b1,d1,c1,b_denom_1) + do j=js,je + do i=is,ie + h_tr = hold(i,j,1) + GV%H_subroundoff + b1(i) = 1.0 / (h_tr + ent(i,j,2)) + d1(i) = h_tr * b1(i) + T(i,j,1) = (b1(i)*h_tr)*T(i,j,1) + S(i,j,1) = (b1(i)*h_tr)*S(i,j,1) + enddo + do k=2,G%ke ; do i=is,ie + c1(i,k) = ent(i,j,K) * b1(i) + h_tr = hold(i,j,k) + GV%H_subroundoff + b_denom_1 = h_tr + d1(i)*ent(i,j,K) + b1(i) = 1.0 / (b_denom_1 + ent(i,j,K+1)) + d1(i) = b_denom_1 * b1(i) + T(i,j,k) = b1(i) * (h_tr*T(i,j,k) + ent(i,j,K)*T(i,j,k-1)) + S(i,j,k) = b1(i) * (h_tr*S(i,j,k) + ent(i,j,K)*S(i,j,k-1)) + enddo ; enddo + do k=G%ke-1,1,-1 ; do i=is,ie + T(i,j,k) = T(i,j,k) + c1(i,k+1)*T(i,j,k+1) + S(i,j,k) = S(i,j,k) + c1(i,k+1)*S(i,j,k+1) + enddo ; enddo + enddo +end subroutine triDiagTS_Eulerian + + !> This subroutine calculates u_h and v_h (velocities at thickness !! points), optionally using the entrainment amounts passed in as arguments. -subroutine find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, ea, eb) +subroutine find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, ea, eb, zero_mix) type(ocean_grid_type), intent(in) :: 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 @@ -449,6 +499,9 @@ subroutine find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, ea, eb) optional, intent(in) :: eb !< The amount of fluid entrained from the layer !! below within this time step [H ~> m or kg m-2]. !! Omitting eb is the same as setting it to 0. + logical, optional, intent(in) :: zero_mix !< If true, do the calculation of u_h and + !! v_h as though ea and eb were being supplied with + !! uniformly zero values. ! local variables real :: b_denom_1 ! The first term in the denominator of b1 [H ~> m or kg m-2]. @@ -459,7 +512,7 @@ subroutine find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, ea, eb) real :: a_e(SZI_(G)), a_w(SZI_(G)) ! velocity points, ~1/2 in the open ! ocean, nondimensional. real :: sum_area, Idenom - logical :: mix_vertically + logical :: mix_vertically, zero_mixing integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke call cpu_clock_begin(id_clock_uv_at_h) @@ -469,6 +522,8 @@ subroutine find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, ea, eb) if (present(ea) .neqv. present(eb)) call MOM_error(FATAL, & "find_uv_at_h: Either both ea and eb or neither one must be present "// & "in call to find_uv_at_h.") + zero_mixing = .false. ; if (present(zero_mix)) zero_mixing = zero_mix + if (zero_mixing) mix_vertically = .false. !$OMP parallel do default(none) shared(is,ie,js,je,G,GV,mix_vertically,h,h_neglect, & !$OMP eb,u_h,u,v_h,v,nz,ea) & !$OMP private(sum_area,Idenom,a_w,a_e,a_s,a_n,b_denom_1,b1,d1,c1) @@ -515,6 +570,17 @@ subroutine find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, ea, eb) u_h(i,j,k) = u_h(i,j,k) + c1(i,k+1)*u_h(i,j,k+1) v_h(i,j,k) = v_h(i,j,k) + c1(i,k+1)*v_h(i,j,k+1) enddo ; enddo + elseif (zero_mixing) then + do i=is,ie + b1(i) = 1.0 / (h(i,j,1) + h_neglect) + u_h(i,j,1) = (h(i,j,1)*b1(i)) * (a_e(i)*u(I,j,1) + a_w(i)*u(I-1,j,1)) + v_h(i,j,1) = (h(i,j,1)*b1(i)) * (a_n(i)*v(i,J,1) + a_s(i)*v(i,J-1,1)) + enddo + do k=2,nz ; do i=is,ie + b1(i) = 1.0 / (h(i,j,k) + h_neglect) + u_h(i,j,k) = (h(i,j,k) * (a_e(i)*u(I,j,k) + a_w(i)*u(I-1,j,k))) * b1(i) + v_h(i,j,k) = (h(i,j,k) * (a_n(i)*v(i,J,k) + a_s(i)*v(i,J-1,k))) * b1(i) + enddo ; enddo else do k=1,nz ; do i=is,ie u_h(i,j,k) = a_e(i)*u(I,j,k) + a_w(i)*u(I-1,j,k) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index e5fe6724c5..ca9e8a45e6 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -13,9 +13,8 @@ module MOM_diabatic_driver use MOM_CVMix_ddiff, only : CVMix_ddiff_is_used use MOM_diabatic_aux, only : diabatic_aux_init, diabatic_aux_end, diabatic_aux_CS use MOM_diabatic_aux, only : make_frazil, adjust_salt, differential_diffuse_T_S, triDiagTS -use MOM_diabatic_aux, only : find_uv_at_h, diagnoseMLDbyDensityDifference, applyBoundaryFluxesInOut -use MOM_diabatic_aux, only : diagnoseMLDbyEnergy -use MOM_diabatic_aux, only : set_pen_shortwave +use MOM_diabatic_aux, only : triDiagTS_Eulerian, find_uv_at_h, diagnoseMLDbyDensityDifference +use MOM_diabatic_aux, only : applyBoundaryFluxesInOut, diagnoseMLDbyEnergy, set_pen_shortwave use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr use MOM_diag_mediator, only : diag_ctrl, time_type, diag_update_remap_grids use MOM_diag_mediator, only : diag_ctrl, query_averaging_enabled, enable_averages, disable_averaging @@ -62,7 +61,7 @@ module MOM_diabatic_driver use MOM_ALE_sponge, only : apply_ALE_sponge, ALE_sponge_CS use MOM_time_manager, only : time_type, real_to_time, operator(-), operator(<=) use MOM_tracer_flow_control, only : call_tracer_column_fns, tracer_flow_control_CS -use MOM_tracer_diabatic, only : tracer_vertdiff +use MOM_tracer_diabatic, only : tracer_vertdiff, tracer_vertdiff_Eulerian use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs, vertvisc_type, accel_diag_ptrs use MOM_variables, only : cont_diag_ptrs, MOM_thermovar_chksum, p3d @@ -573,10 +572,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim if (CS%use_kappa_shear .or. CS%use_CVMix_shear) then if (CS%use_geothermal) then - ! The presence of ea_s causes find_uv_at_h to use a tridiagonal solver, - ! which changes answers at the level of roundoff because ((A*B / A) /= B). - ent_s(:,:,:) = 0.0 - call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, ent_s(:,:,1), ent_s(:,:,1)) + call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, zero_mix=.true.) else call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US) endif @@ -902,10 +898,10 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim ent_t(i,j,K) = ent_s(i,j,K) ; ent_t(i,j,K+1) = ent_s(i,j,K+1) enddo ; enddo ; enddo if (CS%tracer_tridiag) then - call tracer_vertdiff(hold, ent_t(:,:,1:nz), ent_t(:,:,2:nz+1), dt, tv%T, G, GV) - call tracer_vertdiff(hold, ent_s(:,:,1:nz), ent_s(:,:,2:nz+1), dt, tv%S, G, GV) + call tracer_vertdiff_Eulerian(hold, ent_t, dt, tv%T, G, GV) + call tracer_vertdiff_Eulerian(hold, ent_s, dt, tv%S, G, GV) else - call triDiagTS(G, GV, is, ie, js, je, hold, ent_s(:,:,1:nz), ent_s(:,:,2:nz+1), tv%T, tv%S) + call triDiagTS_Eulerian(G, GV, is, ie, js, je, hold, ent_s, tv%T, tv%S) endif ! diagnose temperature, salinity, heat, and salt tendencies @@ -940,12 +936,12 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim if (CS%id_Kd_salt > 0) call post_data(CS%id_Kd_salt, Kd_salt, CS%diag) if (CS%id_Kd_ePBL > 0) call post_data(CS%id_Kd_ePBL, Kd_ePBL, CS%diag) - if (CS%id_ea > 0) call post_data(CS%id_ea, ent_s(:,:,1:nz), CS%diag) - if (CS%id_eb > 0) call post_data(CS%id_eb, ent_s(:,:,2:nz+1), CS%diag) - if (CS%id_ea_t > 0) call post_data(CS%id_ea_t, ent_t(:,:,1:nz), CS%diag) - if (CS%id_eb_t > 0) call post_data(CS%id_eb_t, ent_t(:,:,2:nz+1), CS%diag) - if (CS%id_ea_s > 0) call post_data(CS%id_ea_s, ent_s(:,:,1:nz), CS%diag) - if (CS%id_eb_s > 0) call post_data(CS%id_eb_s, ent_s(:,:,2:nz+1), CS%diag) + if (CS%id_ea > 0) call post_data(CS%id_ea, ent_s(:,:,1:nz), CS%diag) + if (CS%id_eb > 0) call post_data(CS%id_eb, ent_s(:,:,2:nz+1), CS%diag) + if (CS%id_ea_t > 0) call post_data(CS%id_ea_t, ent_t(:,:,1:nz), CS%diag) + if (CS%id_eb_t > 0) call post_data(CS%id_eb_t, ent_t(:,:,2:nz+1), CS%diag) + if (CS%id_ea_s > 0) call post_data(CS%id_ea_s, ent_s(:,:,1:nz), CS%diag) + if (CS%id_eb_s > 0) call post_data(CS%id_eb_s, ent_s(:,:,2:nz+1), CS%diag) Idt = 1.0 / dt if (CS%id_Tdif > 0) then do j=js,je ; do i=is,ie @@ -1183,10 +1179,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, if (CS%use_kappa_shear .or. CS%use_CVMix_shear) then if (CS%use_geothermal) then - ! The presence of ea_s causes find_uv_at_h to use a tridiagonal solver, - ! which changes answers at the level of roundoff because ((A*B / A) /= B). - ! Note that ent_s(:,:,:) = 0.0 - call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, ent_s(:,:,1:nz), ent_s(:,:,1:nz)) + call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, zero_mix=.true.) else call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US) endif @@ -1441,8 +1434,8 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, "Kd_salt (diabatic_ALE)") ! Changes T and S via the tridiagonal solver; no change to h - call tracer_vertdiff(h, ent_t(:,:,1:nz), ent_t(:,:,2:nz+1), dt, tv%T, G, GV) - call tracer_vertdiff(h, ent_s(:,:,1:nz), ent_s(:,:,2:nz+1), dt, tv%S, G, GV) + call tracer_vertdiff_Eulerian(h, ent_t, dt, tv%T, G, GV) + call tracer_vertdiff_Eulerian(h, ent_s, dt, tv%S, G, GV) ! In ALE-mode, layer thicknesses do not change. Therefore, we can use h below if (CS%diabatic_diff_tendency_diag) then diff --git a/src/tracer/MOM_tracer_diabatic.F90 b/src/tracer/MOM_tracer_diabatic.F90 index ec7c025db0..6b9a12f696 100644 --- a/src/tracer/MOM_tracer_diabatic.F90 +++ b/src/tracer/MOM_tracer_diabatic.F90 @@ -13,7 +13,7 @@ module MOM_tracer_diabatic implicit none ; private #include -public tracer_vertdiff +public tracer_vertdiff, tracer_vertdiff_Eulerian public applyTracerBoundaryFluxesInOut contains @@ -41,6 +41,7 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & real, dimension(SZI_(G),SZJ_(G)), optional,intent(in) :: btm_flux !< The (negative upward) bottom flux of the !! tracer in [CU kg m-2 T-1 ~> CU kg m-2 s-1] or !! [CU H ~> CU m or CU kg m-2] if + !! convert_flux_in is .false. real, dimension(SZI_(G),SZJ_(G)), optional,intent(inout) :: btm_reservoir !< amount of tracer in a bottom reservoir !! [CU kg m-2]; formerly [CU m] real, optional,intent(in) :: sink_rate !< rate at which the tracer sinks @@ -70,7 +71,6 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & !! in roundoff and can be neglected [H ~> m or kg m-2]. logical :: convert_flux = .true. - integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -84,20 +84,17 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & h_neglect = GV%H_subroundoff sink_dist = 0.0 if (present(sink_rate)) sink_dist = (dt*sink_rate) * GV%m_to_H -!$OMP parallel default(none) shared(is,ie,js,je,sfc_src,btm_src,sfc_flux,dt,G,GV,btm_flux, & -!$OMP sink_rate,btm_reservoir,nz,sink_dist,ea, & -!$OMP h_old,convert_flux,h_neglect,eb,tr) & -!$OMP private(sink,h_minus_dsink,b_denom_1,b1,d1,h_tr,c1) -!$OMP do + !$OMP parallel default(shared) private(sink,h_minus_dsink,b_denom_1,b1,d1,h_tr,c1) + !$OMP do do j=js,je; do i=is,ie ; sfc_src(i,j) = 0.0 ; btm_src(i,j) = 0.0 ; enddo ; enddo if (present(sfc_flux)) then if (convert_flux) then -!$OMP do + !$OMP do do j = js, je; do i = is,ie sfc_src(i,j) = (sfc_flux(i,j)*dt) * GV%kg_m2_to_H enddo ; enddo else -!$OMP do + !$OMP do do j = js, je; do i = is,ie sfc_src(i,j) = sfc_flux(i,j) enddo ; enddo @@ -105,12 +102,12 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & endif if (present(btm_flux)) then if (convert_flux) then -!$OMP do + !$OMP do do j = js, je; do i = is,ie btm_src(i,j) = (btm_flux(i,j)*dt) * GV%kg_m2_to_H enddo ; enddo else -!$OMP do + !$OMP do do j = js, je; do i = is,ie btm_src(i,j) = btm_flux(i,j) enddo ; enddo @@ -118,7 +115,7 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & endif if (present(sink_rate)) then -!$OMP do + !$OMP do do j=js,je ! Find the sinking rates at all interfaces, limiting them if necesary ! so that the characteristics do not cross within a timestep. @@ -186,7 +183,7 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & endif ; enddo ; enddo enddo else -!$OMP do + !$OMP do do j=js,je do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then h_tr = h_old(i,j,1) + h_neglect @@ -217,11 +214,211 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & endif ; enddo ; enddo enddo endif - -!$OMP end parallel + !$OMP end parallel end subroutine tracer_vertdiff + +!> This subroutine solves a tridiagonal equation for the final tracer concentrations after +!! Eulerian mixing, and possibly sinking or surface and bottom sources, are applied. The sinking +!! is implemented with an fully implicit upwind advection scheme. Alternate time units can be +!! used for the timestep, surface and bottom fluxes and sink_rate provided they are all consistent. +subroutine tracer_vertdiff_Eulerian(h_old, ent, dt, tr, G, GV, & + sfc_flux, btm_flux, btm_reservoir, sink_rate, convert_flux_in) + type(ocean_grid_type), intent(in) :: G !< ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_old !< layer thickness before entrainment + !! [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: ent !< Amount of fluid mixed across interfaces + !! [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: tr !< tracer concentration in concentration units [CU] + real, intent(in) :: dt !< amount of time covered by this call [T ~> s] + real, dimension(SZI_(G),SZJ_(G)), optional,intent(in) :: sfc_flux !< surface flux of the tracer in units of + !! [CU kg m-2 T-1 ~> CU kg m-2 s-1] or + !! [CU H ~> CU m or CU kg m-2] if + !! convert_flux_in is .false. + real, dimension(SZI_(G),SZJ_(G)), optional,intent(in) :: btm_flux !< The (negative upward) bottom flux of the + !! tracer in [CU kg m-2 T-1 ~> CU kg m-2 s-1] or + !! [CU H ~> CU m or CU kg m-2] if + !! convert_flux_in is .false. + real, dimension(SZI_(G),SZJ_(G)), optional,intent(inout) :: btm_reservoir !< amount of tracer in a bottom reservoir + !! [CU kg m-2]; formerly [CU m] + real, optional,intent(in) :: sink_rate !< rate at which the tracer sinks + !! [m T-1 ~> m s-1] + logical, optional,intent(in) :: convert_flux_in !< True if the specified sfc_flux needs + !! to be integrated in time + + ! local variables + real :: sink_dist !< The distance the tracer sinks in a time step [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G)) :: & + sfc_src, & !< The time-integrated surface source of the tracer [CU H ~> CU m or CU kg m-2]. + btm_src !< The time-integrated bottom source of the tracer [CU H ~> CU m or CU kg m-2]. + real, dimension(SZI_(G)) :: & + b1, & !< b1 is used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1]. + d1 !! d1=1-c1 is used by the tridiagonal solver, nondimensional. + real :: c1(SZI_(G),SZK_(GV)) !< c1 is used by the tridiagonal solver [nondim]. + real :: h_minus_dsink(SZI_(G),SZK_(GV)) !< The layer thickness minus the + !! difference in sinking rates across the layer [H ~> m or kg m-2]. + !! By construction, 0 <= h_minus_dsink < h_work. + real :: sink(SZI_(G),SZK_(GV)+1) !< The tracer's sinking distances at the + !! interfaces, limited to prevent characteristics from + !! crossing within a single timestep [H ~> m or kg m-2]. + real :: b_denom_1 !< The first term in the denominator of b1 [H ~> m or kg m-2]. + real :: h_tr !< h_tr is h at tracer points with a h_neglect added to + !! ensure positive definiteness [H ~> m or kg m-2]. + real :: h_neglect !< A thickness that is so small it is usually lost + !! in roundoff and can be neglected [H ~> m or kg m-2]. + logical :: convert_flux = .true. + + + integer :: i, j, k, is, ie, js, je, nz + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + + if (nz == 1) then + call MOM_error(WARNING, "MOM_tracer_diabatic.F90, tracer_vertdiff called "//& + "with only one vertical level") + return + endif + + if (present(convert_flux_in)) convert_flux = convert_flux_in + h_neglect = GV%H_subroundoff + sink_dist = 0.0 + if (present(sink_rate)) sink_dist = (dt*sink_rate) * GV%m_to_H + !$OMP parallel default(shared) private(sink,h_minus_dsink,b_denom_1,b1,d1,h_tr,c1) + !$OMP do + do j=js,je; do i=is,ie ; sfc_src(i,j) = 0.0 ; btm_src(i,j) = 0.0 ; enddo ; enddo + if (present(sfc_flux)) then + if (convert_flux) then + !$OMP do + do j = js, je; do i = is,ie + sfc_src(i,j) = (sfc_flux(i,j)*dt) * GV%kg_m2_to_H + enddo ; enddo + else + !$OMP do + do j = js, je; do i = is,ie + sfc_src(i,j) = sfc_flux(i,j) + enddo ; enddo + endif + endif + if (present(btm_flux)) then + if (convert_flux) then + !$OMP do + do j = js, je; do i = is,ie + btm_src(i,j) = (btm_flux(i,j)*dt) * GV%kg_m2_to_H + enddo ; enddo + else + !$OMP do + do j = js, je; do i = is,ie + btm_src(i,j) = btm_flux(i,j) + enddo ; enddo + endif + endif + + if (present(sink_rate)) then + !$OMP do + do j=js,je + ! Find the sinking rates at all interfaces, limiting them if necesary + ! so that the characteristics do not cross within a timestep. + ! If a non-constant sinking rate were used, that would be incorprated + ! here. + if (present(btm_reservoir)) then + do i=is,ie ; sink(i,nz+1) = sink_dist ; enddo + do k=2,nz ; do i=is,ie + sink(i,K) = sink_dist ; h_minus_dsink(i,k) = h_old(i,j,k) + enddo ; enddo + else + do i=is,ie ; sink(i,nz+1) = 0.0 ; enddo + ! Find the limited sinking distance at the interfaces. + do k=nz,2,-1 ; do i=is,ie + if (sink(i,K+1) >= sink_dist) then + sink(i,K) = sink_dist + h_minus_dsink(i,k) = h_old(i,j,k) + (sink(i,K+1) - sink(i,K)) + elseif (sink(i,K+1) + h_old(i,j,k) < sink_dist) then + sink(i,K) = sink(i,K+1) + h_old(i,j,k) + h_minus_dsink(i,k) = 0.0 + else + sink(i,K) = sink_dist + h_minus_dsink(i,k) = (h_old(i,j,k) + sink(i,K+1)) - sink(i,K) + endif + enddo ; enddo + endif + do i=is,ie + sink(i,1) = 0.0 ; h_minus_dsink(i,1) = (h_old(i,j,1) + sink(i,2)) + enddo + + ! Now solve the tridiagonal equation for the tracer concentrations. + do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then + b_denom_1 = h_minus_dsink(i,1) + ent(i,j,1) + h_neglect + b1(i) = 1.0 / (b_denom_1 + ent(i,j,2)) + d1(i) = b_denom_1 * b1(i) + h_tr = h_old(i,j,1) + h_neglect + tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + sfc_src(i,j) + endif ; enddo + do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then + c1(i,k) = ent(i,j,K) * b1(i) + b_denom_1 = h_minus_dsink(i,k) + d1(i) * (ent(i,j,K) + sink(i,K)) + & + h_neglect + b1(i) = 1.0 / (b_denom_1 + ent(i,j,K+1)) + d1(i) = b_denom_1 * b1(i) + h_tr = h_old(i,j,k) + h_neglect + tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + & + (ent(i,j,K) + sink(i,K)) * tr(i,j,k-1)) + endif ; enddo ; enddo + do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then + c1(i,nz) = ent(i,j,nz) * b1(i) + b_denom_1 = h_minus_dsink(i,nz) + d1(i) * (ent(i,j,nz) + sink(i,nz)) + & + h_neglect + b1(i) = 1.0 / (b_denom_1 + ent(i,j,nz+1)) + h_tr = h_old(i,j,nz) + h_neglect + tr(i,j,nz) = b1(i) * ((h_tr * tr(i,j,nz) + btm_src(i,j)) + & + (ent(i,j,nz) + sink(i,nz)) * tr(i,j,nz-1)) + endif ; enddo + if (present(btm_reservoir)) then ; do i=is,ie ; if (G%mask2dT(i,j)>0.5) then + btm_reservoir(i,j) = btm_reservoir(i,j) + & + (sink(i,nz+1)*tr(i,j,nz)) * GV%H_to_kg_m2 + endif ; enddo ; endif + + do k=nz-1,1,-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then + tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1) + endif ; enddo ; enddo + enddo + else + !$OMP do + do j=js,je + do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + h_tr = h_old(i,j,1) + h_neglect + b_denom_1 = h_tr + ent(i,j,1) + b1(i) = 1.0 / (b_denom_1 + ent(i,j,2)) + d1(i) = h_tr * b1(i) + tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + sfc_src(i,j) + endif + enddo + do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + c1(i,k) = ent(i,j,K) * b1(i) + h_tr = h_old(i,j,k) + h_neglect + b_denom_1 = h_tr + d1(i) * ent(i,j,K) + b1(i) = 1.0 / (b_denom_1 + ent(i,j,K+1)) + d1(i) = b_denom_1 * b1(i) + tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + ent(i,j,K) * tr(i,j,k-1)) + endif ; enddo ; enddo + do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + c1(i,nz) = ent(i,j,nz) * b1(i) + h_tr = h_old(i,j,nz) + h_neglect + b_denom_1 = h_tr + d1(i)*ent(i,j,nz) + b1(i) = 1.0 / ( b_denom_1 + ent(i,j,nz+1)) + tr(i,j,nz) = b1(i) * (( h_tr * tr(i,j,nz) + btm_src(i,j)) + & + ent(i,j,nz) * tr(i,j,nz-1)) + endif ; enddo + do k=nz-1,1,-1 ; do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1) + endif ; enddo ; enddo + enddo + endif + !$OMP end parallel + +end subroutine tracer_vertdiff_Eulerian + + !> This routine is modeled after applyBoundaryFluxesInOut in MOM_diabatic_aux.F90 !! NOTE: Please note that in this routine sfc_flux gets set to zero to ensure that the surface !! flux of the tracer does not get applied again during a subsequent call to tracer_vertdif From 15cf6e0ec655f172fda63129267348ae7747db74 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 25 Nov 2020 11:42:52 -0500 Subject: [PATCH 03/12] Refactoring and code cleanup in MOM_ALE_sponge.F90 Code refactoring and cleanup in MOM_ALE_sponge.F90 to eliminate array-syntax calculations, remove unnecessary variables, avoid duplicated large-stride memory copies, remove nonsensical comments and fix indents. All answers are bitwise identical. --- .../vertical/MOM_ALE_sponge.F90 | 204 ++++++++---------- src/parameterizations/vertical/MOM_sponge.F90 | 2 - 2 files changed, 94 insertions(+), 112 deletions(-) diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index edd35d61f2..7ef0877321 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -162,8 +162,6 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ #include "version_variable.h" character(len=40) :: mdl = "MOM_sponge" ! This module's name. logical :: use_sponge - real, allocatable, dimension(:,:,:) :: data_hu !< thickness at u points [H ~> m or kg m-2] - real, allocatable, dimension(:,:,:) :: data_hv !< thickness at v points [H ~> m or kg m-2] real, allocatable, dimension(:,:) :: Iresttime_u !< inverse of the restoring time at u points [T-1 ~> s-1] real, allocatable, dimension(:,:) :: Iresttime_v !< inverse of the restoring time at v points [T-1 ~> s-1] logical :: bndExtrapolation = .true. ! If true, extrapolate boundaries @@ -259,43 +257,40 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ "The total number of columns where sponges are applied at h points.", like_default=.true.) if (CS%sponge_uv) then - - allocate(data_hu(G%isdB:G%iedB,G%jsd:G%jed,nz_data)) ; data_hu(:,:,:) = 0.0 - allocate(data_hv(G%isd:G%ied,G%jsdB:G%jedB,nz_data)) ; data_hv(:,:,:) = 0.0 allocate(Iresttime_u(G%isdB:G%iedB,G%jsd:G%jed)) ; Iresttime_u(:,:) = 0.0 allocate(Iresttime_v(G%isd:G%ied,G%jsdB:G%jedB)) ; Iresttime_v(:,:) = 0.0 ! u points CS%num_col_u = 0 ; !CS%fldno_u = 0 do j=CS%jsc,CS%jec ; do I=CS%iscB,CS%iecB - data_hu(I,j,:) = 0.5 * (data_h(i,j,:) + data_h(i+1,j,:)) - Iresttime_u(I,j) = 0.5 * (Iresttime(i,j) + Iresttime(i+1,j)) - if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) & - CS%num_col_u = CS%num_col_u + 1 + Iresttime_u(I,j) = 0.5 * (Iresttime(i,j) + Iresttime(i+1,j)) + if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) CS%num_col_u = CS%num_col_u + 1 enddo ; enddo if (CS%num_col_u > 0) then - allocate(CS%Iresttime_col_u(CS%num_col_u)) ; CS%Iresttime_col_u(:) = 0.0 - allocate(CS%col_i_u(CS%num_col_u)) ; CS%col_i_u(:) = 0 - allocate(CS%col_j_u(CS%num_col_u)) ; CS%col_j_u(:) = 0 - - ! pass indices, restoring time to the CS structure - col = 1 - do j=CS%jsc,CS%jec ; do I=CS%iscB,CS%iecB - if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) then - CS%col_i_u(col) = i ; CS%col_j_u(col) = j - CS%Iresttime_col_u(col) = Iresttime_u(i,j) - col = col + 1 - endif - enddo ; enddo - - ! same for total number of arbritary layers and correspondent data - - allocate(CS%Ref_hu%p(CS%nz_data,CS%num_col_u)) - do col=1,CS%num_col_u ; do K=1,CS%nz_data - CS%Ref_hu%p(K,col) = data_hu(CS%col_i_u(col),CS%col_j_u(col),K) - enddo ; enddo + allocate(CS%Iresttime_col_u(CS%num_col_u)) ; CS%Iresttime_col_u(:) = 0.0 + allocate(CS%col_i_u(CS%num_col_u)) ; CS%col_i_u(:) = 0 + allocate(CS%col_j_u(CS%num_col_u)) ; CS%col_j_u(:) = 0 + + ! Store the column indices and restoring rates in the CS structure + col = 1 + do j=CS%jsc,CS%jec ; do I=CS%iscB,CS%iecB + if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) then + CS%col_i_u(col) = I ; CS%col_j_u(col) = j + CS%Iresttime_col_u(col) = Iresttime_u(I,j) + col = col + 1 + endif + enddo ; enddo + + ! same for total number of arbritary layers and correspondent data + allocate(CS%Ref_hu%p(CS%nz_data,CS%num_col_u)) + do col=1,CS%num_col_u + I = CS%col_i_u(col) ; j = CS%col_j_u(col) + do k=1,CS%nz_data + CS%Ref_hu%p(k,col) = 0.5 * (data_h(i,j,k) + data_h(i+1,j,k)) + enddo + enddo endif total_sponge_cols_u = CS%num_col_u call sum_across_PEs(total_sponge_cols_u) @@ -305,10 +300,8 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ ! v points CS%num_col_v = 0 ; !CS%fldno_v = 0 do J=CS%jscB,CS%jecB; do i=CS%isc,CS%iec - data_hv(i,J,:) = 0.5 * (data_h(i,j,:) + data_h(i,j+1,:)) Iresttime_v(i,J) = 0.5 * (Iresttime(i,j) + Iresttime(i,j+1)) - if ((Iresttime_v(i,J)>0.0) .and. (G%mask2dCv(i,J)>0)) & - CS%num_col_v = CS%num_col_v + 1 + if ((Iresttime_v(i,J)>0.0) .and. (G%mask2dCv(i,J)>0)) CS%num_col_v = CS%num_col_v + 1 enddo ; enddo if (CS%num_col_v > 0) then @@ -329,9 +322,12 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ ! same for total number of arbritary layers and correspondent data allocate(CS%Ref_hv%p(CS%nz_data,CS%num_col_v)) - do col=1,CS%num_col_v ; do K=1,CS%nz_data - CS%Ref_hv%p(K,col) = data_hv(CS%col_i_v(col),CS%col_j_v(col),K) - enddo ; enddo + do col=1,CS%num_col_v + i = CS%col_i_v(col) ; J = CS%col_j_v(col) + do k=1,CS%nz_data + CS%Ref_hv%p(k,col) = 0.5 * (data_h(i,j,k) + data_h(i,j+1,k)) + enddo + enddo endif total_sponge_cols_v = CS%num_col_v call sum_across_PEs(total_sponge_cols_v) @@ -473,7 +469,7 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, param_file, CS) if ((Iresttime(i,j)>0.0) .and. (G%mask2dT(i,j)>0)) then CS%col_i(col) = i ; CS%col_j(col) = j CS%Iresttime_col(col) = Iresttime(i,j) - col = col +1 + col = col + 1 endif enddo ; enddo endif @@ -505,7 +501,7 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, param_file, CS) if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) then CS%col_i_u(col) = i ; CS%col_j_u(col) = j CS%Iresttime_col_u(col) = Iresttime_u(i,j) - col = col +1 + col = col + 1 endif enddo ; enddo ! same for total number of arbritary layers and correspondent data @@ -531,7 +527,7 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, param_file, CS) if ((Iresttime_v(i,J)>0.0) .and. (G%mask2dCv(i,J)>0)) then CS%col_i_v(col) = i ; CS%col_j_v(col) = j CS%Iresttime_col_v(col) = Iresttime_v(i,j) - col = col +1 + col = col + 1 endif enddo ; enddo endif @@ -800,8 +796,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) real :: m_to_Z ! A unit conversion factor from m to Z. real, allocatable, dimension(:) :: tmp_val2 ! data values on the original grid real, dimension(SZK_(G)) :: tmp_val1 ! data values remapped to model grid - real :: hu(SZIB_(G), SZJ_(G), SZK_(G)) ! A temporary array for h at u pts - real :: hv(SZI_(G), SZJB_(G), SZK_(G)) ! A temporary array for h at v pts + real, dimension(SZK_(G)) :: h_col ! A column of thicknesses at h, u or v points [H ~> m or kg m-2] real, allocatable, dimension(:,:,:) :: sp_val ! A temporary array for fields real, allocatable, dimension(:,:,:) :: mask_z ! A temporary array for field mask at h pts real, dimension(:), allocatable :: hsrc ! Source thicknesses [Z ~> m]. @@ -830,47 +825,45 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) if (.not. present(Time)) & call MOM_error(FATAL,"apply_ALE_sponge: No time information provided") do m=1,CS%fldno - nz_data = CS%Ref_val(m)%nz_data - allocate(sp_val(G%isd:G%ied,G%jsd:G%jed,1:nz_data)) - allocate(mask_z(G%isd:G%ied,G%jsd:G%jed,1:nz_data)) - sp_val(:,:,:) = 0.0 - mask_z(:,:,:) = 0.0 - call horiz_interp_and_extrap_tracer(CS%Ref_val(m)%id, Time, 1.0, G, sp_val, mask_z, z_in, & + nz_data = CS%Ref_val(m)%nz_data + allocate(sp_val(G%isd:G%ied,G%jsd:G%jed,1:nz_data)) ; sp_val(:,:,:) = 0.0 + allocate(mask_z(G%isd:G%ied,G%jsd:G%jed,1:nz_data)) ; mask_z(:,:,:) = 0.0 + call horiz_interp_and_extrap_tracer(CS%Ref_val(m)%id, Time, 1.0, G, sp_val, mask_z, z_in, & z_edges_in, missing_value, .true., .false., .false., & spongeOnGrid=CS%SpongeDataOngrid, m_to_Z=US%m_to_Z, & answers_2018=CS%hor_regrid_answers_2018) - allocate( hsrc(nz_data) ) - allocate( tmpT1d(nz_data) ) - do c=1,CS%num_col - i = CS%col_i(c) ; j = CS%col_j(c) - CS%Ref_val(m)%p(1:nz_data,c) = sp_val(i,j,1:nz_data) - ! Build the source grid - zTopOfCell = 0. ; zBottomOfCell = 0. ; nPoints = 0; hsrc(:) = 0.0; tmpT1d(:) = -99.9 - do k=1,nz_data - if (mask_z(CS%col_i(c),CS%col_j(c),k) == 1.0) then - zBottomOfCell = -min( z_edges_in(k+1), G%bathyT(CS%col_i(c),CS%col_j(c)) ) - tmpT1d(k) = sp_val(CS%col_i(c),CS%col_j(c),k) - elseif (k>1) then - zBottomOfCell = -G%bathyT(CS%col_i(c),CS%col_j(c)) - tmpT1d(k) = tmpT1d(k-1) - else ! This next block should only ever be reached over land - tmpT1d(k) = -99.9 - endif - hsrc(k) = zTopOfCell - zBottomOfCell - if (hsrc(k)>0.) nPoints = nPoints + 1 - zTopOfCell = zBottomOfCell ! Bottom becomes top for next value of k - enddo - ! In case data is deeper than model - hsrc(nz_data) = hsrc(nz_data) + ( zTopOfCell + G%bathyT(CS%col_i(c),CS%col_j(c)) ) - CS%Ref_val(m)%h(1:nz_data,c) = GV%Z_to_H*hsrc(1:nz_data) - CS%Ref_val(m)%p(1:nz_data,c) = tmpT1d(1:nz_data) - do k=2,nz_data - ! if (mask_z(i,j,k)==0.) & - if (CS%Ref_val(m)%h(k,c) <= 0.001*GV%m_to_H) & - ! some confusion here about why the masks are not correct returning from horiz_interp - ! reverting to using a minimum thickness criteria - CS%Ref_val(m)%p(k,c) = CS%Ref_val(m)%p(k-1,c) - enddo + allocate( hsrc(nz_data) ) + allocate( tmpT1d(nz_data) ) + do c=1,CS%num_col + i = CS%col_i(c) ; j = CS%col_j(c) + do k=1,nz_data ; CS%Ref_val(m)%p(k,c) = sp_val(i,j,k) ; enddo + ! Build the source grid + zTopOfCell = 0. ; zBottomOfCell = 0. ; nPoints = 0 ; hsrc(:) = 0. ; tmpT1d(:) = -99.9 + do k=1,nz_data + if (mask_z(CS%col_i(c),CS%col_j(c),k) == 1.0) then + zBottomOfCell = -min( z_edges_in(k+1), G%bathyT(CS%col_i(c),CS%col_j(c)) ) + tmpT1d(k) = sp_val(CS%col_i(c),CS%col_j(c),k) + elseif (k>1) then + zBottomOfCell = -G%bathyT(CS%col_i(c),CS%col_j(c)) + tmpT1d(k) = tmpT1d(k-1) + else ! This next block should only ever be reached over land + tmpT1d(k) = -99.9 + endif + hsrc(k) = zTopOfCell - zBottomOfCell + if (hsrc(k)>0.) nPoints = nPoints + 1 + zTopOfCell = zBottomOfCell ! Bottom becomes top for next value of k + enddo + ! In case data is deeper than model + hsrc(nz_data) = hsrc(nz_data) + ( zTopOfCell + G%bathyT(CS%col_i(c),CS%col_j(c)) ) + CS%Ref_val(m)%h(1:nz_data,c) = GV%Z_to_H*hsrc(1:nz_data) + CS%Ref_val(m)%p(1:nz_data,c) = tmpT1d(1:nz_data) + do k=2,nz_data + ! if (mask_z(i,j,k)==0.) & + if (CS%Ref_val(m)%h(k,c) <= 0.001*GV%m_to_H) & + ! some confusion here about why the masks are not correct returning from horiz_interp + ! reverting to using a minimum thickness criteria + CS%Ref_val(m)%p(k,c) = CS%Ref_val(m)%p(k-1,c) + enddo enddo deallocate(sp_val, mask_z, hsrc, tmpT1d) enddo @@ -879,23 +872,22 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) endif allocate(tmp_val2(nz_data)) - do m=1,CS%fldno - do c=1,CS%num_col -! c is an index for the next 3 lines but a multiplier for the rest of the loop -! Therefore we use c as per C code and increment the index where necessary. - i = CS%col_i(c) ; j = CS%col_j(c) - damp = dt * CS%Iresttime_col(c) - I1pdamp = 1.0 / (1.0 + damp) + do c=1,CS%num_col + i = CS%col_i(c) ; j = CS%col_j(c) + damp = dt * CS%Iresttime_col(c) + I1pdamp = 1.0 / (1.0 + damp) + do k=1,nz ; h_col(k) = h(i,j,k) ; enddo + do m=1,CS%fldno tmp_val2(1:nz_data) = CS%Ref_val(m)%p(1:nz_data,c) if (CS%time_varying_sponges) then call remapping_core_h(CS%remap_cs, nz_data, CS%Ref_val(m)%h(1:nz_data,c), tmp_val2, & - CS%nz, h(i,j,:), tmp_val1, h_neglect, h_neglect_edge) + CS%nz, h_col, tmp_val1, h_neglect, h_neglect_edge) else - call remapping_core_h(CS%remap_cs,nz_data, CS%Ref_h%p(1:nz_data,c), tmp_val2, & - CS%nz, h(i,j,:), tmp_val1, h_neglect, h_neglect_edge) + call remapping_core_h(CS%remap_cs, nz_data, CS%Ref_h%p(1:nz_data,c), tmp_val2, & + CS%nz, h_col, tmp_val1, h_neglect, h_neglect_edge) endif - !Backward Euler method - CS%var(m)%p(i,j,1:CS%nz) = I1pdamp * (CS%var(m)%p(i,j,1:CS%nz) + tmp_val1 * damp) + ! Backward Euler method + do k=1,CS%nz ; CS%var(m)%p(i,j,k) = I1pdamp * (CS%var(m)%p(i,j,k) + tmp_val1(k) * damp) ; enddo enddo enddo @@ -907,10 +899,6 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) !enddo if (CS%sponge_uv) then - ! u points - do j=CS%jsc,CS%jec; do I=CS%iscB,CS%iecB; do k=1,nz - hu(I,j,k) = 0.5 * (h(i,j,k) + h(i+1,j,k)) - enddo ; enddo ; enddo if (CS%time_varying_sponges) then if (.not. present(Time)) & call MOM_error(FATAL,"apply_ALE_sponge: No time information provided") @@ -925,10 +913,8 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) ! call pass_var(sp_val,G%Domain) ! call pass_var(mask_z,G%Domain) do c=1,CS%num_col - ! c is an index for the next 3 lines but a multiplier for the rest of the loop - ! Therefore we use c as per C code and increment the index where necessary. i = CS%col_i(c) ; j = CS%col_j(c) - CS%Ref_val_u%p(1:nz_data,c) = sp_val(i,j,1:nz_data) + do k=1,nz_data ; CS%Ref_val_u%p(k,c) = sp_val(i,j,k) ; enddo enddo deallocate (sp_val, mask_z) @@ -945,10 +931,8 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) ! call pass_var(mask_z,G%Domain) do c=1,CS%num_col - ! c is an index for the next 3 lines but a multiplier for the rest of the loop - ! Therefore we use c as per C code and increment the index where necessary. i = CS%col_i(c) ; j = CS%col_j(c) - CS%Ref_val_v%p(1:nz_data,c) = sp_val(i,j,1:nz_data) + do k=1,nz_data ; CS%Ref_val_v%p(k,c) = sp_val(i,j,k) ; enddo enddo deallocate (sp_val, mask_z) @@ -957,42 +941,42 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) nz_data = CS%nz_data endif + ! u points do c=1,CS%num_col_u - i = CS%col_i_u(c) ; j = CS%col_j_u(c) + I = CS%col_i_u(c) ; j = CS%col_j_u(c) damp = dt * CS%Iresttime_col_u(c) I1pdamp = 1.0 / (1.0 + damp) + do k=1,nz ; h_col(k) = 0.5 * (h(i,j,k) + h(i+1,j,k)) ; enddo if (CS%time_varying_sponges) nz_data = CS%Ref_val(m)%nz_data tmp_val2(1:nz_data) = CS%Ref_val_u%p(1:nz_data,c) if (CS%time_varying_sponges) then call remapping_core_h(CS%remap_cs, nz_data, CS%Ref_val_u%h(:,c), tmp_val2, & - CS%nz, hu(i,j,:), tmp_val1, h_neglect, h_neglect_edge) + CS%nz, h_col, tmp_val1, h_neglect, h_neglect_edge) else call remapping_core_h(CS%remap_cs, nz_data, CS%Ref_hu%p(:,c), tmp_val2, & - CS%nz, hu(i,j,:), tmp_val1, h_neglect, h_neglect_edge) + CS%nz, h_col, tmp_val1, h_neglect, h_neglect_edge) endif !Backward Euler method - CS%var_u%p(i,j,:) = I1pdamp * (CS%var_u%p(i,j,:) + tmp_val1 * damp) + do k=1,CS%nz ; CS%var_u%p(I,j,k) = I1pdamp * (CS%var_u%p(I,j,k) + tmp_val1(k) * damp) ; enddo enddo ! v points - do J=CS%jscB,CS%jecB; do i=CS%isc,CS%iec; do k=1,nz - hv(i,J,k) = 0.5 * (h(i,j,k) + h(i,j+1,k)) - enddo ; enddo ; enddo do c=1,CS%num_col_v i = CS%col_i_v(c) ; j = CS%col_j_v(c) damp = dt * CS%Iresttime_col_v(c) I1pdamp = 1.0 / (1.0 + damp) + do k=1,nz ; h_col(k) = 0.5 * (h(i,j,k) + h(i,j+1,k)) ; enddo tmp_val2(1:nz_data) = CS%Ref_val_v%p(1:nz_data,c) if (CS%time_varying_sponges) then call remapping_core_h(CS%remap_cs, CS%nz_data, CS%Ref_val_v%h(:,c), tmp_val2, & - CS%nz, hv(i,j,:), tmp_val1, h_neglect, h_neglect_edge) + CS%nz, h_col, tmp_val1, h_neglect, h_neglect_edge) else call remapping_core_h(CS%remap_cs, CS%nz_data, CS%Ref_hv%p(:,c), tmp_val2, & - CS%nz, hv(i,j,:), tmp_val1, h_neglect, h_neglect_edge) + CS%nz, h_col, tmp_val1, h_neglect, h_neglect_edge) endif !Backward Euler method - CS%var_v%p(i,j,:) = I1pdamp * (CS%var_v%p(i,j,:) + tmp_val1 * damp) + do k=1,CS%nz ; CS%var_u%p(i,J,k) = I1pdamp * (CS%var_u%p(i,J,k) + tmp_val1(k) * damp) ; enddo enddo endif diff --git a/src/parameterizations/vertical/MOM_sponge.F90 b/src/parameterizations/vertical/MOM_sponge.F90 index 4566abcef7..dcd0ac4e02 100644 --- a/src/parameterizations/vertical/MOM_sponge.F90 +++ b/src/parameterizations/vertical/MOM_sponge.F90 @@ -477,8 +477,6 @@ subroutine apply_sponge(h, dt, G, GV, US, ea, eb, CS, Rcv_ml) endif do c=1,CS%num_col -! c is an index for the next 3 lines but a multiplier for the rest of the loop -! Therefore we use c as per C code and increment the index where necessary. i = CS%col_i(c) ; j = CS%col_j(c) damp = dt * CS%Iresttime_col(c) From 8be13be646e419e53a537801f7e9cbce0ce51145 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 25 Nov 2020 11:43:11 -0500 Subject: [PATCH 04/12] Resize internal args in MOM_bulk_mixed_layer.F90 Changed the vertical extent of several arguments to internal subroutines in MOM_bulk_mixed_layer.F90 to avoid possibly forcing array copies and simplify the calls. All answers are bitwise identical. --- .../vertical/MOM_bulk_mixed_layer.F90 | 93 +++++++++---------- 1 file changed, 42 insertions(+), 51 deletions(-) diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index 079655f787..3bea0d9937 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -482,11 +482,10 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C if (CS%ML_resort) then if (id_clock_resort>0) call cpu_clock_begin(id_clock_resort) if (CS%ML_presort_nz_conv_adj > 0) & - call convective_adjustment(h(:,1:), u, v, R0(:,1:), Rcv(:,1:), T(:,1:), & - S(:,1:), eps, d_eb, dKE_CA, cTKE, j, G, GV, US, CS, & - CS%ML_presort_nz_conv_adj) + call convective_adjustment(h, u, v, R0, Rcv, T, S, eps, d_eb, dKE_CA, cTKE, j, G, GV, & + US, CS, CS%ML_presort_nz_conv_adj) - call sort_ML(h(:,1:), R0(:,1:), eps, G, GV, CS, ksort) + call sort_ML(h, R0, eps, G, GV, CS, ksort) if (id_clock_resort>0) call cpu_clock_end(id_clock_resort) else do k=1,nz ; do i=is,ie ; ksort(i,k) = k ; enddo ; enddo @@ -495,8 +494,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C ! Undergo instantaneous entrainment into the buffer layers and mixed layers ! to remove hydrostatic instabilities. Any water that is lighter than ! currently in the mixed or buffer layer is entrained. - call convective_adjustment(h(:,1:), u, v, R0(:,1:), Rcv(:,1:), T(:,1:), & - S(:,1:), eps, d_eb, dKE_CA, cTKE, j, G, GV, US, CS) + call convective_adjustment(h, u, v, R0, Rcv, T, S, eps, d_eb, dKE_CA, cTKE, j, G, GV, US, CS) do i=is,ie ; h_CA(i) = h(i,1) ; enddo if (id_clock_adjustment>0) call cpu_clock_end(id_clock_adjustment) @@ -535,18 +533,15 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C ! Pen_SW_bnd = components to penetrative shortwave radiation call extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & CS%H_limit_fluxes, CS%use_river_heat_content, CS%use_calving_heat_content, & - h(:,1:), T(:,1:), netMassInOut, netMassOut, Net_heat, Net_salt, Pen_SW_bnd,& + h(:,1:), T(:,1:), netMassInOut, netMassOut, Net_heat, Net_salt, Pen_SW_bnd, & tv, aggregate_FW_forcing) ! This subroutine causes the mixed layer to entrain to depth of free convection. - call mixedlayer_convection(h(:,1:), d_eb, htot, Ttot, Stot, uhtot, vhtot, & - R0_tot, Rcv_tot, u, v, T(:,1:), S(:,1:), & - R0(:,1:), Rcv(:,1:), eps, & - dR0_dT, dRcv_dT, dR0_dS, dRcv_dS, & - netMassInOut, netMassOut, Net_heat, Net_salt, & - nsw, Pen_SW_bnd, opacity_band, Conv_En, & - dKE_FC, j, ksort, G, GV, US, CS, tv, fluxes, dt, & - aggregate_FW_forcing) + call mixedlayer_convection(h, d_eb, htot, Ttot, Stot, uhtot, vhtot, R0_tot, Rcv_tot, & + u, v, T, S, R0, Rcv, eps, dR0_dT, dRcv_dT, dR0_dS, dRcv_dS, & + netMassInOut, netMassOut, Net_heat, Net_salt, & + nsw, Pen_SW_bnd, opacity_band, Conv_En, dKE_FC, & + j, ksort, G, GV, US, CS, tv, fluxes, dt, aggregate_FW_forcing) if (id_clock_conv>0) call cpu_clock_end(id_clock_conv) @@ -563,8 +558,8 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C j, ksort, G, GV, US, CS) ! Here the mechanically driven entrainment occurs. - call mechanical_entrainment(h(:,1:), d_eb, htot, Ttot, Stot, uhtot, vhtot, & - R0_tot, Rcv_tot, u, v, T(:,1:), S(:,1:), R0(:,1:), Rcv(:,1:), eps, dR0_dT, dRcv_dT, & + call mechanical_entrainment(h, d_eb, htot, Ttot, Stot, uhtot, vhtot, & + R0_tot, Rcv_tot, u, v, T, S, R0, Rcv, eps, dR0_dT, dRcv_dT, & cMKE, Idt_diag, nsw, Pen_SW_bnd, opacity_band, TKE, & Idecay_len_TKE, j, ksort, G, GV, US, CS) @@ -803,43 +798,39 @@ end subroutine bulkmixedlayer !! is lighter than currently in the mixed- or buffer- layer is entrained. subroutine convective_adjustment(h, u, v, R0, Rcv, T, S, eps, d_eb, & dKE_CA, cTKE, j, G, GV, US, CS, nz_conv) - 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(inout) :: h !< Layer thickness [H ~> m or kg m-2]. + 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),SZK0_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]. !! The units of h are referred to as H below. - real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: u !< Zonal velocities interpolated to h + real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: u !< Zonal velocities interpolated to h !! points [L T-1 ~> m s-1]. - real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: v !< Zonal velocities interpolated to h + real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: v !< Zonal velocities interpolated to h !! points [L T-1 ~> m s-1]. - real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: T !< Layer temperatures [degC]. - real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: S !< Layer salinities [ppt]. - real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: R0 !< Potential density referenced to + real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: R0 !< Potential density referenced to !! surface pressure [R ~> kg m-3]. - real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: Rcv !< The coordinate defining potential + real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: Rcv !< The coordinate defining potential !! density [R ~> kg m-3]. - real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: d_eb !< The downward increase across a layer + real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: T !< Layer temperatures [degC]. + real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: S !< Layer salinities [ppt]. + real, dimension(SZI_(G),SZK_(GV)), intent(in) :: eps !< The negligibly small amount of water + !! that will be left in each layer [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: d_eb !< The downward increase across a layer !! in the entrainment from below [H ~> m or kg m-2]. !! Positive values go with mass gain by !! a layer. - real, dimension(SZI_(G),SZK_(GV)), intent(in) :: eps !< The negligibly small amount of water - !! that will be left in each layer [H ~> m or kg m-2]. - real, dimension(SZI_(G),SZK_(GV)), intent(out) :: dKE_CA !< The vertically integrated change in + real, dimension(SZI_(G),SZK_(GV)), intent(out) :: dKE_CA !< The vertically integrated change in !! kinetic energy due to convective !! adjustment [Z L2 T-2 ~> m3 s-2]. - real, dimension(SZI_(G),SZK_(GV)), intent(out) :: cTKE !< The buoyant turbulent kinetic energy + real, dimension(SZI_(G),SZK_(GV)), intent(out) :: cTKE !< The buoyant turbulent kinetic energy !! source due to convective adjustment !! [Z L2 T-2 ~> m3 s-2]. - integer, intent(in) :: j !< The j-index to work on. - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - type(bulkmixedlayer_CS), pointer :: CS !< The control structure for this module. - integer, optional, intent(in) :: nz_conv !< If present, the number of layers + integer, intent(in) :: j !< The j-index to work on. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(bulkmixedlayer_CS), pointer :: CS !< The control structure for this module. + integer, optional, intent(in) :: nz_conv !< If present, the number of layers !! over which to do convective adjustment !! (perhaps CS%nkml). -! This subroutine does instantaneous convective entrainment into the buffer -! layers and mixed layers to remove hydrostatic instabilities. Any water that -! is lighter than currently in the mixed- or buffer- layer is entrained. - ! Local variables real, dimension(SZI_(G)) :: & htot, & ! The total depth of the layers being considered for @@ -942,7 +933,7 @@ subroutine mixedlayer_convection(h, d_eb, htot, Ttot, Stot, uhtot, vhtot, & aggregate_FW_forcing) 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)), & + real, dimension(SZI_(G),SZK0_(GV)), & intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]. !! The units of h are referred to as H below. real, dimension(SZI_(G),SZK_(GV)), & @@ -966,14 +957,14 @@ subroutine mixedlayer_convection(h, d_eb, htot, Ttot, Stot, uhtot, vhtot, & intent(in) :: u !< Zonal velocities interpolated to h points [L T-1 ~> m s-1]. real, dimension(SZI_(G),SZK_(GV)), & intent(in) :: v !< Zonal velocities interpolated to h points [L T-1 ~> m s-1]. - real, dimension(SZI_(G),SZK_(GV)), & + real, dimension(SZI_(G),SZK0_(GV)), & intent(in) :: T !< Layer temperatures [degC]. - real, dimension(SZI_(G),SZK_(GV)), & + real, dimension(SZI_(G),SZK0_(GV)), & intent(in) :: S !< Layer salinities [ppt]. - real, dimension(SZI_(G),SZK_(GV)), & + real, dimension(SZI_(G),SZK0_(GV)), & intent(in) :: R0 !< Potential density referenced to !! surface pressure [R ~> kg m-3]. - real, dimension(SZI_(G),SZK_(GV)), & + real, dimension(SZI_(G),SZK0_(GV)), & intent(in) :: Rcv !< The coordinate defining potential !! density [R ~> kg m-3]. real, dimension(SZI_(G),SZK_(GV)), & @@ -1503,7 +1494,7 @@ subroutine mechanical_entrainment(h, d_eb, htot, Ttot, Stot, uhtot, vhtot, & type(ocean_grid_type), intent(in) :: 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 - real, dimension(SZI_(G),SZK_(GV)), & + real, dimension(SZI_(G),SZK0_(GV)), & intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]. real, dimension(SZI_(G),SZK_(GV)), & intent(inout) :: d_eb !< The downward increase across a layer in the @@ -1526,14 +1517,14 @@ subroutine mechanical_entrainment(h, d_eb, htot, Ttot, Stot, uhtot, vhtot, & intent(in) :: u !< Zonal velocities interpolated to h points [L T-1 ~> m s-1]. real, dimension(SZI_(G),SZK_(GV)), & intent(in) :: v !< Zonal velocities interpolated to h points [L T-1 ~> m s-1]. - real, dimension(SZI_(G),SZK_(GV)), & + real, dimension(SZI_(G),SZK0_(GV)), & intent(in) :: T !< Layer temperatures [degC]. - real, dimension(SZI_(G),SZK_(GV)), & + real, dimension(SZI_(G),SZK0_(GV)), & intent(in) :: S !< Layer salinities [ppt]. - real, dimension(SZI_(G),SZK_(GV)), & + real, dimension(SZI_(G),SZK0_(GV)), & intent(in) :: R0 !< Potential density referenced to !! surface pressure [R ~> kg m-3]. - real, dimension(SZI_(G),SZK_(GV)), & + real, dimension(SZI_(G),SZK0_(GV)), & intent(in) :: Rcv !< The coordinate defining potential !! density [R ~> kg m-3]. real, dimension(SZI_(G),SZK_(GV)), & @@ -1840,8 +1831,8 @@ end subroutine mechanical_entrainment subroutine sort_ML(h, R0, eps, G, GV, CS, ksort) 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) :: h !< Layer thickness [H ~> m or kg m-2]. - real, dimension(SZI_(G),SZK_(GV)), intent(in) :: R0 !< The potential density used to sort + real, dimension(SZI_(G),SZK0_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZK0_(GV)), intent(in) :: R0 !< The potential density used to sort !! the layers [R ~> kg m-3]. real, dimension(SZI_(G),SZK_(GV)), intent(in) :: eps !< The (small) thickness that must !! remain in each layer [H ~> m or kg m-2]. From 3ceb13dfc67982c77f7d9a7b6ea76223c7ffd654 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 25 Nov 2020 11:43:38 -0500 Subject: [PATCH 05/12] +Eliminate diabatic diagnostics that are always 0 Eliminated register_diag_field and post_data calls for several diagnostics in modes where they are always all zeros. This will change some output fields and the list of available diags in some cases, so it is expected to fail the automated diagnostic regression tests. All solutions are bitwise identical. --- .../vertical/MOM_diabatic_driver.F90 | 144 +++++++----------- 1 file changed, 55 insertions(+), 89 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index ca9e8a45e6..2fb6a27542 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -173,14 +173,17 @@ module MOM_diabatic_driver !>@{ Diagnostic IDs integer :: id_cg1 = -1 ! diag handle for mode-1 speed integer, allocatable, dimension(:) :: id_cn ! diag handle for all mode speeds - integer :: id_wd = -1, id_ea = -1, id_eb = -1 ! used by layer diabatic - integer :: id_dudt_dia = -1, id_dvdt_dia = -1, id_ea_s = -1, id_eb_s = -1 - integer :: id_ea_t = -1, id_eb_t = -1 - integer :: id_Kd_heat = -1, id_Kd_salt = -1, id_Kd_interface = -1, id_Kd_ePBL = -1 - integer :: id_Tdif = -1, id_Tadv = -1, id_Sdif = -1, id_Sadv = -1 - integer :: id_MLD_003 = -1, id_MLD_0125 = -1, id_MLD_user = -1, id_mlotstsq = -1 - integer :: id_MLD_EN1 = -1, id_MLD_EN2 = -1, id_MLD_EN3= -1 - integer :: id_subMLN2 = -1 + integer :: id_ea = -1, id_eb = -1 ! used by layer diabatic + integer :: id_ea_t = -1, id_eb_t = -1, id_ea_s = -1, id_eb_s = -1 + integer :: id_Kd_heat = -1, id_Kd_salt = -1, id_Kd_int = -1, id_Kd_ePBL = -1 + integer :: id_Tdif = -1, id_Sdif = -1, id_Tadv = -1, id_Sadv = -1 + ! These are handles to diagnostics related to the mixed layer properties. + integer :: id_MLD_003 = -1, id_MLD_0125 = -1, id_MLD_user = -1, id_mlotstsq = -1 + integer :: id_MLD_EN1 = -1, id_MLD_EN2 = -1, id_MLD_EN3 = -1, id_subMLN2 = -1 + + ! These are handles to diatgnostics that are only available in non-ALE layered mode. + integer :: id_wd = -1 + integer :: id_dudt_dia = -1, id_dvdt_dia = -1 integer :: id_hf_dudt_dia_2d = -1, id_hf_dvdt_dia_2d = -1 ! diagnostic for fields prior to applying diapycnal physics @@ -501,16 +504,10 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim ! Kd_int [Z2 T-1 ~> m2 s-1]. Kd_extra_S , & ! The extra diffusivity of salinity due to double diffusion relative to ! Kd_int [Z2 T-1 ~> m2 s-1]. - zeros_h, & ! An array of zeros to handle diagnostics that should be removed. Kd_ePBL, & ! test array of diapycnal diffusivities at interfaces [Z2 T-1 ~> m2 s-1] Tdif_flx, & ! diffusive diapycnal heat flux across interfaces [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1] Sdif_flx ! diffusive diapycnal salt flux across interfaces [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1] - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: & - zeros_u ! An array of zeros for u-point diagnostics that should be removed. - real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: & - zeros_v ! An array of zeros for v-point diagnostics that should be removed. - logical :: in_boundary(SZI_(G)) ! True if there are no massive layers below, ! where massive is defined as sufficiently thick that ! the no-flux boundary conditions have not restricted @@ -931,10 +928,10 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim call diag_update_remap_grids(CS%diag) ! Diagnose the diapycnal diffusivities and other related quantities. - if (CS%id_Kd_interface > 0) call post_data(CS%id_Kd_interface, Kd_int, CS%diag) - if (CS%id_Kd_heat > 0) call post_data(CS%id_Kd_heat, Kd_heat, CS%diag) - if (CS%id_Kd_salt > 0) call post_data(CS%id_Kd_salt, Kd_salt, CS%diag) - if (CS%id_Kd_ePBL > 0) call post_data(CS%id_Kd_ePBL, Kd_ePBL, CS%diag) + if (CS%id_Kd_int > 0) call post_data(CS%id_Kd_int, Kd_int, CS%diag) + if (CS%id_Kd_heat > 0) call post_data(CS%id_Kd_heat, Kd_heat, CS%diag) + if (CS%id_Kd_salt > 0) call post_data(CS%id_Kd_salt, Kd_salt, CS%diag) + if (CS%id_Kd_ePBL > 0) call post_data(CS%id_Kd_ePBL, Kd_ePBL, CS%diag) if (CS%id_ea > 0) call post_data(CS%id_ea, ent_s(:,:,1:nz), CS%diag) if (CS%id_eb > 0) call post_data(CS%id_eb, ent_s(:,:,2:nz+1), CS%diag) @@ -1039,18 +1036,6 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim endif endif ! CS%use_sponge - !! Diagnostics that are always zero and should be removed from this routine. - if ((CS%id_wd > 0) .or. (CS%id_Tadv > 0) .or. (CS%id_Sadv > 0)) zeros_h(:,:,:) = 0.0 - if ((CS%id_dudt_dia > 0) .or. (CS%id_hf_dudt_dia_2d > 0)) zeros_u(:,:,:) = 0.0 - if ((CS%id_dvdt_dia > 0) .or. (CS%id_hf_dvdt_dia_2d > 0)) zeros_v(:,:,:) = 0.0 - if (CS%id_wd > 0) call post_data(CS%id_wd, zeros_h, CS%diag) - if (CS%id_Tadv > 0) call post_data(CS%id_Tadv, zeros_h, CS%diag) - if (CS%id_Sadv > 0) call post_data(CS%id_Sadv, zeros_h, CS%diag) - if (CS%id_dudt_dia > 0) call post_data(CS%id_dudt_dia, zeros_u, CS%diag) - if (CS%id_dvdt_dia > 0) call post_data(CS%id_dvdt_dia, zeros_v, CS%diag) - if (CS%id_hf_dudt_dia_2d > 0) call post_data(CS%id_hf_dudt_dia_2d, zeros_u(:,:,1), CS%diag) - if (CS%id_hf_dvdt_dia_2d > 0) call post_data(CS%id_hf_dvdt_dia_2d, zeros_v(:,:,1), CS%diag) - call disable_averaging(CS%diag) if (showCallTree) call callTree_leave("diabatic_ALE_legacy()") @@ -1110,13 +1095,8 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, Kd_extra_S , & ! The extra diffusivity of salinity due to double diffusion relative to ! Kd_int returned from set_diffusivity [Z2 T-1 ~> m2 s-1]. Kd_ePBL, & ! boundary layer or convective diapycnal diffusivities at interfaces [Z2 T-1 ~> m2 s-1] - zeros_h, & ! An array of zeros for h-point diagnostics that should be removed. Tdif_flx, & ! diffusive diapycnal heat flux across interfaces [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1] Sdif_flx ! diffusive diapycnal salt flux across interfaces [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1] - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: & - zeros_u ! An array of zeros for u-point diagnostics that should be removed. - real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: & - zeros_v ! An array of zeros for v-point diagnostics that should be removed. logical :: in_boundary(SZI_(G)) ! True if there are no massive layers below, ! where massive is defined as sufficiently thick that @@ -1209,7 +1189,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, endif ! Store the diagnosed typical diffusivity at interfaces. - if (CS%id_Kd_interface > 0) call post_data(CS%id_Kd_interface, Kd_heat, CS%diag) + if (CS%id_Kd_int > 0) call post_data(CS%id_Kd_int, Kd_heat, CS%diag) ! Set diffusivities for heat and salt separately, and possibly change the meaning of Kd_heat. if (CS%double_diffuse) then @@ -1548,17 +1528,6 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, call pass_var(visc%Kv_slow, G%Domain, To_All+Omit_Corners, halo=1) call cpu_clock_end(id_clock_pass) - !! Diagnostics that are always zero and should be removed from this routine. - if ((CS%id_dudt_dia > 0) .or. (CS%id_hf_dudt_dia_2d > 0)) zeros_u(:,:,:) = 0.0 - if ((CS%id_dvdt_dia > 0) .or. (CS%id_hf_dvdt_dia_2d > 0)) zeros_v(:,:,:) = 0.0 - if ((CS%id_Tadv > 0) .or. (CS%id_Sadv > 0)) zeros_h(:,:,:) = 0.0 - if (CS%id_dudt_dia > 0) call post_data(CS%id_dudt_dia, zeros_u, CS%diag) - if (CS%id_dvdt_dia > 0) call post_data(CS%id_dvdt_dia, zeros_v, CS%diag) - if (CS%id_hf_dudt_dia_2d > 0) call post_data(CS%id_hf_dudt_dia_2d, zeros_u(:,:,1), CS%diag) - if (CS%id_hf_dvdt_dia_2d > 0) call post_data(CS%id_hf_dvdt_dia_2d, zeros_v(:,:,1), CS%diag) - if (CS%id_Tadv > 0) call post_data(CS%id_Tadv, zeros_h, CS%diag) - if (CS%id_Sadv > 0) call post_data(CS%id_Sadv, zeros_h, CS%diag) - call disable_averaging(CS%diag) if (showCallTree) call callTree_leave("diabatic_ALE()") @@ -1627,7 +1596,6 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e ! Kd_int [Z2 T-1 ~> m2 s-1]. Kd_extra_S , & ! The extra diffusivity of salinity due to double diffusion relative to ! Kd_int [Z2 T-1 ~> m2 s-1]. - zeros_h, & ! An array of zeros to handle diagnostics that should be removed. Tdif_flx, & ! diffusive diapycnal heat flux across interfaces [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1] Tadv_flx, & ! advective diapycnal heat flux across interfaces [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1] Sdif_flx, & ! diffusive diapycnal salt flux across interfaces [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1] @@ -2479,14 +2447,12 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e endif ! Diagnose the diapycnal diffusivities and other related quantities. - if (CS%id_Kd_interface > 0) call post_data(CS%id_Kd_interface, Kd_int, CS%diag) - if (CS%id_Kd_heat > 0) call post_data(CS%id_Kd_heat, Kd_heat, CS%diag) - if (CS%id_Kd_salt > 0) call post_data(CS%id_Kd_salt, Kd_salt, CS%diag) - if (CS%id_Kd_ePBL > 0) zeros_h(:,:,:) = 0.0 - if (CS%id_Kd_ePBL > 0) call post_data(CS%id_Kd_ePBL, zeros_h, CS%diag) + if (CS%id_Kd_int > 0) call post_data(CS%id_Kd_int, Kd_int, CS%diag) + if (CS%id_Kd_heat > 0) call post_data(CS%id_Kd_heat, Kd_heat, CS%diag) + if (CS%id_Kd_salt > 0) call post_data(CS%id_Kd_salt, Kd_salt, CS%diag) - if (CS%id_ea > 0) call post_data(CS%id_ea, ea, CS%diag) - if (CS%id_eb > 0) call post_data(CS%id_eb, eb, CS%diag) + if (CS%id_ea > 0) call post_data(CS%id_ea, ea, CS%diag) + if (CS%id_eb > 0) call post_data(CS%id_eb, eb, CS%diag) if (CS%id_dudt_dia > 0) call post_data(CS%id_dudt_dia, ADp%du_dt_dia, CS%diag) if (CS%id_dvdt_dia > 0) call post_data(CS%id_dvdt_dia, ADp%dv_dt_dia, CS%diag) @@ -3051,32 +3017,31 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di 'Layer entrainment from above per timestep', 'm', conversion=GV%H_to_m) CS%id_eb = register_diag_field('ocean_model', 'eb', diag%axesTL, Time, & 'Layer entrainment from below per timestep', 'm', conversion=GV%H_to_m) - !### if (.not.CS%useALEalgorithm) then - CS%id_wd = register_diag_field('ocean_model', 'wd', diag%axesTi, Time, & - 'Diapycnal velocity', 'm s-1', conversion=GV%H_to_m) - if (CS%id_wd > 0) call safe_alloc_ptr(CDp%diapyc_vel,isd,ied,jsd,jed,nz+1) - - CS%id_dudt_dia = register_diag_field('ocean_model', 'dudt_dia', diag%axesCuL, Time, & - 'Zonal Acceleration from Diapycnal Mixing', 'm s-2', conversion=US%L_T2_to_m_s2) - CS%id_dvdt_dia = register_diag_field('ocean_model', 'dvdt_dia', diag%axesCvL, Time, & - 'Meridional Acceleration from Diapycnal Mixing', 'm s-2', conversion=US%L_T2_to_m_s2) - - CS%id_hf_dudt_dia_2d = register_diag_field('ocean_model', 'hf_dudt_dia_2d', diag%axesCu1, Time, & - 'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Diapycnal Mixing', & - 'm s-2', conversion=US%L_T2_to_m_s2) - if (CS%id_hf_dudt_dia_2d > 0) then - call safe_alloc_ptr(ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz) - call safe_alloc_ptr(ADp%du_dt_dia,IsdB,IedB,jsd,jed,nz) - endif + if (.not.CS%useALEalgorithm) then + CS%id_wd = register_diag_field('ocean_model', 'wd', diag%axesTi, Time, & + 'Diapycnal velocity', 'm s-1', conversion=GV%H_to_m) + if (CS%id_wd > 0) call safe_alloc_ptr(CDp%diapyc_vel,isd,ied,jsd,jed,nz+1) - CS%id_hf_dvdt_dia_2d = register_diag_field('ocean_model', 'hf_dvdt_dia_2d', diag%axesCv1, Time, & - 'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Diapycnal Mixing', & - 'm s-2', conversion=US%L_T2_to_m_s2) - if (CS%id_hf_dvdt_dia_2d > 0) then - call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz) - call safe_alloc_ptr(ADp%dv_dt_dia,isd,ied,JsdB,JedB,nz) + CS%id_dudt_dia = register_diag_field('ocean_model', 'dudt_dia', diag%axesCuL, Time, & + 'Zonal Acceleration from Diapycnal Mixing', 'm s-2', conversion=US%L_T2_to_m_s2) + CS%id_dvdt_dia = register_diag_field('ocean_model', 'dvdt_dia', diag%axesCvL, Time, & + 'Meridional Acceleration from Diapycnal Mixing', 'm s-2', conversion=US%L_T2_to_m_s2) + + CS%id_hf_dudt_dia_2d = register_diag_field('ocean_model', 'hf_dudt_dia_2d', diag%axesCu1, Time, & + 'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Diapycnal Mixing', & + 'm s-2', conversion=US%L_T2_to_m_s2) + if (CS%id_hf_dudt_dia_2d > 0) call safe_alloc_ptr(ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz) + + CS%id_hf_dvdt_dia_2d = register_diag_field('ocean_model', 'hf_dvdt_dia_2d', diag%axesCv1, Time, & + 'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Diapycnal Mixing', & + 'm s-2', conversion=US%L_T2_to_m_s2) + if (CS%id_hf_dvdt_dia_2d > 0) call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + + if ((CS%id_dudt_dia > 0) .or. (CS%id_hf_dudt_dia_2d > 0)) & + call safe_alloc_ptr(ADp%du_dt_dia,IsdB,IedB,jsd,jed,nz) + if ((CS%id_dvdt_dia > 0) .or. (CS%id_hf_dvdt_dia_2d > 0)) & + call safe_alloc_ptr(ADp%dv_dt_dia,isd,ied,JsdB,JedB,nz) endif - !### endif if (CS%use_int_tides) then CS%id_cg1 = register_diag_field('ocean_model', 'cn1', diag%axesT1, & @@ -3095,15 +3060,19 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di CS%id_Tdif = register_diag_field('ocean_model',"Tflx_dia_diff", diag%axesTi, & Time, "Diffusive diapycnal temperature flux across interfaces", & "degC m s-1", conversion=GV%H_to_m*US%s_to_T) - CS%id_Tadv = register_diag_field('ocean_model',"Tflx_dia_adv", diag%axesTi, & - Time, "Advective diapycnal temperature flux across interfaces", & - "degC m s-1", conversion=GV%H_to_m*US%s_to_T) + if (.not.CS%useALEalgorithm) then + CS%id_Tadv = register_diag_field('ocean_model',"Tflx_dia_adv", diag%axesTi, & + Time, "Advective diapycnal temperature flux across interfaces", & + "degC m s-1", conversion=GV%H_to_m*US%s_to_T) + endif CS%id_Sdif = register_diag_field('ocean_model',"Sflx_dia_diff", diag%axesTi, & Time, "Diffusive diapycnal salnity flux across interfaces", & "psu m s-1", conversion=GV%H_to_m*US%s_to_T) - CS%id_Sadv = register_diag_field('ocean_model',"Sflx_dia_adv", diag%axesTi, & - Time, "Advective diapycnal salnity flux across interfaces", & - "psu m s-1", conversion=GV%H_to_m*US%s_to_T) + if (.not.CS%useALEalgorithm) then + CS%id_Sadv = register_diag_field('ocean_model',"Sflx_dia_adv", diag%axesTi, & + Time, "Advective diapycnal salnity flux across interfaces", & + "psu m s-1", conversion=GV%H_to_m*US%s_to_T) + endif CS%id_MLD_003 = register_diag_field('ocean_model', 'MLD_003', diag%axesT1, Time, & 'Mixed layer depth (delta rho = 0.03)', 'm', conversion=US%Z_to_m, & cmor_field_name='mlotst', cmor_long_name='Ocean Mixed Layer Thickness Defined by Sigma T', & @@ -3151,9 +3120,6 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di "stratification at the base of the mixed layer.", & units='m', default=50.0, scale=US%m_to_Z) - if (CS%id_dudt_dia > 0) call safe_alloc_ptr(ADp%du_dt_dia,IsdB,IedB,jsd,jed,nz) - if (CS%id_dvdt_dia > 0) call safe_alloc_ptr(ADp%dv_dt_dia,isd,ied,JsdB,JedB,nz) - ! diagnostics for values prior to diabatic and prior to ALE CS%id_u_predia = register_diag_field('ocean_model', 'u_predia', diag%axesCuL, Time, & 'Zonal velocity before diabatic forcing', 'm s-1', conversion=US%L_T_to_m_s) @@ -3173,7 +3139,7 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di !call set_diffusivity_init(Time, G, param_file, diag, CS%set_diff_CSp, CS%int_tide_CSp) - CS%id_Kd_interface = register_diag_field('ocean_model', 'Kd_interface', diag%axesTi, Time, & + 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 CS%id_Kd_ePBL = register_diag_field('ocean_model', 'Kd_ePBL', diag%axesTi, Time, & From bac594cbdf4c09af6460e2d4d9a685587673086c Mon Sep 17 00:00:00 2001 From: Matthew Thompson Date: Wed, 25 Nov 2020 14:39:00 -0500 Subject: [PATCH 06/12] Use flush function --- config_src/solo_driver/MOM_driver.F90 | 2 +- src/diagnostics/MOM_PointAccel.F90 | 4 ++-- src/framework/MOM_domains.F90 | 2 +- src/framework/MOM_write_cputime.F90 | 4 ++-- 4 files changed, 6 insertions(+), 6 deletions(-) diff --git a/config_src/solo_driver/MOM_driver.F90 b/config_src/solo_driver/MOM_driver.F90 index ba52d9c02a..da517b6e9d 100644 --- a/config_src/solo_driver/MOM_driver.F90 +++ b/config_src/solo_driver/MOM_driver.F90 @@ -257,7 +257,7 @@ program MOM_main !$ call omp_set_num_threads(ocean_nthreads) !$OMP PARALLEL !$ write(6,*) "ocean_solo OMPthreading ", fms_affinity_get(), omp_get_thread_num(), omp_get_num_threads() -!$ call flush(6) +!$ flush(6) !$OMP END PARALLEL ! Read ocean_solo restart, which can override settings from the namelist. diff --git a/src/diagnostics/MOM_PointAccel.F90 b/src/diagnostics/MOM_PointAccel.F90 index f6326b06fa..091f88fac2 100644 --- a/src/diagnostics/MOM_PointAccel.F90 +++ b/src/diagnostics/MOM_PointAccel.F90 @@ -390,7 +390,7 @@ subroutine write_u_accel(I, j, um, hin, ADp, CDp, dt_in_T, G, GV, US, CS, vel_rp write(file,'(2/)') - call flush(file) + flush(file) endif end subroutine write_u_accel @@ -722,7 +722,7 @@ subroutine write_v_accel(i, J, vm, hin, ADp, CDp, dt_in_T, G, GV, US, CS, vel_rp write(file,'(2/)') - call flush(file) + flush(file) endif end subroutine write_v_accel diff --git a/src/framework/MOM_domains.F90 b/src/framework/MOM_domains.F90 index 67989f37da..46cc9c526a 100644 --- a/src/framework/MOM_domains.F90 +++ b/src/framework/MOM_domains.F90 @@ -1309,7 +1309,7 @@ subroutine MOM_domains_init(MOM_dom, param_file, symmetric, static_memory, & !$ call fms_affinity_set('OCEAN', ocean_omp_hyper_thread, ocean_nthreads) !$ call omp_set_num_threads(ocean_nthreads) !$ write(6,*) "MOM_domains_mod OMPthreading ", fms_affinity_get(), omp_get_thread_num(), omp_get_num_threads() -!$ call flush(6) +!$ flush(6) !$ endif #endif call log_param(param_file, mdl, "!SYMMETRIC_MEMORY_", MOM_dom%symmetric, & diff --git a/src/framework/MOM_write_cputime.F90 b/src/framework/MOM_write_cputime.F90 index 2c1cb3378a..4ef00707fe 100644 --- a/src/framework/MOM_write_cputime.F90 +++ b/src/framework/MOM_write_cputime.F90 @@ -110,7 +110,7 @@ subroutine MOM_write_cputime_end(CS) ! Flush and close the output files. if (is_root_pe() .and. CS%fileCPU_ascii > 0) then - call flush(CS%fileCPU_ascii) + flush(CS%fileCPU_ascii) call close_file(CS%fileCPU_ascii) endif @@ -200,7 +200,7 @@ subroutine write_cputime(day, n, CS, nmax, call_end) reday, n, (CS%cputime2 / real(CLOCKS_PER_SEC)), & d_cputime / real(CLOCKS_PER_SEC) - call flush(CS%fileCPU_ascii) + flush(CS%fileCPU_ascii) endif CS%previous_calls = CS%previous_calls + 1 From a62867a27913f285d74a5a4f86429b87d36a4f70 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 25 Nov 2020 17:59:55 -0500 Subject: [PATCH 07/12] Allocate a diagnostic pointer to avoid a seg-fault Added a diagnostic pointer allocation statement depending on whether KE_dia is registered to avoid a possible segmentation fault, and corrected an openMP directive, so that the automated test cases pass, apart from expected changes diagnostic outputs. All answers are bitwise identical. --- src/diagnostics/MOM_diagnostics.F90 | 1 + src/parameterizations/vertical/MOM_diabatic_aux.F90 | 6 +++--- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 3936a788d0..1f55801064 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -2225,6 +2225,7 @@ subroutine set_dependent_diagnostics(MIS, ADp, CDp, G, CS) if (associated(CS%KE_dia)) then call safe_alloc_ptr(ADp%du_dt_dia,IsdB,IedB,jsd,jed,nz) call safe_alloc_ptr(ADp%dv_dt_dia,isd,ied,JsdB,JedB,nz) + call safe_alloc_ptr(CDp%diapyc_vel,isd,ied,jsd,jed,nz+1) endif if (associated(CS%uhGM_Rlay)) call safe_alloc_ptr(CDp%uhGM,IsdB,IedB,jsd,jed,nz) diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 40e7f16e76..1683e21fbe 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -524,9 +524,9 @@ subroutine find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, ea, eb, zero_mix) "in call to find_uv_at_h.") zero_mixing = .false. ; if (present(zero_mix)) zero_mixing = zero_mix if (zero_mixing) mix_vertically = .false. -!$OMP parallel do default(none) shared(is,ie,js,je,G,GV,mix_vertically,h,h_neglect, & -!$OMP eb,u_h,u,v_h,v,nz,ea) & -!$OMP private(sum_area,Idenom,a_w,a_e,a_s,a_n,b_denom_1,b1,d1,c1) + !$OMP parallel do default(none) shared(is,ie,js,je,G,GV,mix_vertically,zero_mixing,h, & + !$OMP h_neglect,ea,eb,u_h,u,v_h,v,nz) & + !$OMP private(sum_area,Idenom,a_w,a_e,a_s,a_n,b_denom_1,b1,d1,c1) do j=js,je do i=is,ie sum_area = G%areaCu(I-1,j) + G%areaCu(I,j) From bf420736d97856ed03f6428ab6d6f5a94d724278 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 27 Nov 2020 09:26:01 -0500 Subject: [PATCH 08/12] (*)Use hvel_shelf in vertvisc_coef with ice shelves Explicitly set hvel_shelf in all cases that have shelf points, to avoid using uninitialized arrays when setting CS%h_u and CS%h_v in some cases, and also to simplify the code by reducing the number of calls to find_coupling_coef. All answers in MOM6-examples test cases are bitwise identical, and it seems probable that if there are cases where answers would change that these might fail due to the use of uninitialized variables. --- .../vertical/MOM_vert_friction.F90 | 22 ++++++++----------- 1 file changed, 9 insertions(+), 13 deletions(-) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 7786bf5b46..81f8b29a63 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -822,9 +822,7 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC) enddo if (do_any_shelf) then if (CS%harmonic_visc) then - call find_coupling_coef(a_shelf, hvel, do_i_shelf, h_harm, bbl_thick, & - kv_bbl, z_i, h_ml, dt, j, G, GV, US, CS, & - visc, forces, work_on_u=.true., OBC=OBC, shelf=.true.) + do k=1,nz ; do I=Isq,Ieq ; hvel_shelf(I,k) = hvel(I,k) ; enddo ; enddo else ! Find upwind-biased thickness near the surface. ! Perhaps this needs to be done more carefully, via find_eta. do I=Isq,Ieq ; if (do_i_shelf(I)) then @@ -850,10 +848,10 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC) endif endif ; enddo enddo - call find_coupling_coef(a_shelf, hvel_shelf, do_i_shelf, h_harm, & - bbl_thick, kv_bbl, z_i, h_ml, dt, j, G, GV, US, CS, & - visc, forces, work_on_u=.true., OBC=OBC, shelf=.true.) endif + call find_coupling_coef(a_shelf, hvel_shelf, do_i_shelf, h_harm, bbl_thick, & + kv_bbl, z_i, h_ml, dt, j, G, GV, US, CS, visc, forces, & + work_on_u=.true., OBC=OBC, shelf=.true.) do I=Isq,Ieq ; if (do_i_shelf(I)) CS%a1_shelf_u(I,j) = a_shelf(I,1) ; enddo endif endif @@ -979,7 +977,7 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC) call find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, h_ml, & dt, j, G, GV, US, CS, visc, forces, work_on_u=.false., OBC=OBC) if ( allocated(hML_v)) then - do i=is,ie ; if (do_i(i)) then ; hML_v(i,J) = h_ml(i) ; endif ; enddo + do i=is,ie ; if (do_i(i)) then ; hML_v(i,J) = h_ml(i) ; endif ; enddo endif do_any_shelf = .false. if (associated(forces%frac_shelf_v)) then @@ -990,9 +988,7 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC) enddo if (do_any_shelf) then if (CS%harmonic_visc) then - call find_coupling_coef(a_shelf, hvel, do_i_shelf, h_harm, bbl_thick, & - kv_bbl, z_i, h_ml, dt, j, G, GV, US, CS, visc, & - forces, work_on_u=.false., OBC=OBC, shelf=.true.) + do k=1,nz ; do i=is,ie ; hvel_shelf(i,k) = hvel(i,k) ; enddo ; enddo else ! Find upwind-biased thickness near the surface. ! Perhaps this needs to be done more carefully, via find_eta. do i=is,ie ; if (do_i_shelf(i)) then @@ -1018,10 +1014,10 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC) endif endif ; enddo enddo - call find_coupling_coef(a_shelf, hvel_shelf, do_i_shelf, h_harm, & - bbl_thick, kv_bbl, z_i, h_ml, dt, j, G, GV, US, CS, & - visc, forces, work_on_u=.false., OBC=OBC, shelf=.true.) endif + call find_coupling_coef(a_shelf, hvel_shelf, do_i_shelf, h_harm, bbl_thick, & + kv_bbl, z_i, h_ml, dt, j, G, GV, US, CS, visc, forces, & + work_on_u=.false., OBC=OBC, shelf=.true.) do i=is,ie ; if (do_i_shelf(i)) CS%a1_shelf_v(i,J) = a_shelf(i,1) ; enddo endif endif From 443e7d80395bad339b8c24c3ea8f7765f5853551 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 29 Nov 2020 13:07:02 -0500 Subject: [PATCH 09/12] +Add option to write shorter available diags files Added the option to write out available diags files that avoid having a huge number of redundant entries related to alternative diagnostic grids or module names and derived quantities like spatial averages. In this case, these alternative variants are documented in the available diags files with a brief 'variants' line. A new runtime parameter, FULL_AVAILABLE_DIAGS_LIST, was added to enable this new behavior; by default the available diags files look like they did previously, but this might be the exception to the rule where code changes should not change output by default. Also cleaned up some non-standard code style (like white-space) in MOM_diag_mediator. All answers are bitwise identical, but there is a new entry in the MOM_parameter_doc.all files. --- src/framework/MOM_diag_mediator.F90 | 642 +++++++++++++++------------- 1 file changed, 344 insertions(+), 298 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 28c4c867d7..a929613fbc 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -242,7 +242,10 @@ module MOM_diag_mediator !! This file is open if available_diag_doc_unit is > 0. integer :: chksum_iounit = -1 !< The unit number of a diagnostic documentation file. !! This file is open if available_diag_doc_unit is > 0. - logical :: diag_as_chksum !< If true, log chksums in a text file instead of posting diagnostics + logical :: full_diag_list = .true. !< If true, write an entry for every variant of a + !! diagnostic to the available_diags file; otherwise only + !! write an entry for the primary diagnostics. + logical :: diag_as_chksum !< If true, log chksums in a text file instead of posting diagnostics logical :: grid_space_axes !< If true, diagnostic horizontal coordinates axes are in grid space. ! The following fields are used for the output of the data. integer :: is !< The start i-index of cell centers within the computational domain @@ -606,7 +609,7 @@ subroutine set_axes_info_dsamp(G, GV, param_file, diag_cs, id_zl_native, id_zi_n id_zl = id_zl_native ; id_zi = id_zi_native !Axes group for native downsampled diagnostics do dl=2,MAX_DSAMP_LEV - if(dl .ne. 2) call MOM_error(FATAL, "set_axes_info_dsamp: Downsample level other than 2 is not supported yet!") + if (dl /= 2) call MOM_error(FATAL, "set_axes_info_dsamp: Downsample level other than 2 is not supported yet!") if (G%symmetric) then allocate(gridLonB_dsamp(diag_cs%dsamp(dl)%isgB:diag_cs%dsamp(dl)%iegB)) allocate(gridLatB_dsamp(diag_cs%dsamp(dl)%jsgB:diag_cs%dsamp(dl)%jegB)) @@ -871,7 +874,7 @@ subroutine set_masks_for_axes_dsamp(G, diag_cs) !The downsampled mask is needed for sending out the diagnostics output via diag_manager !The non-downsampled mask is needed for downsampling the diagnostics field do dl=2,MAX_DSAMP_LEV - if(dl .ne. 2) call MOM_error(FATAL, "set_masks_for_axes_dsamp: Downsample level other than 2 is not supported!") + if (dl /= 2) call MOM_error(FATAL, "set_masks_for_axes_dsamp: Downsample level other than 2 is not supported!") do c=1, diag_cs%num_diag_coords ! Level/layer h-points in diagnostic coordinate axes => diag_cs%remap_axesTL(c) @@ -1281,7 +1284,7 @@ subroutine post_data_0d(diag_field_id, field, diag_cs, is_static) if (diag_cs%diag_as_chksum) then call chksum0(locfield, diag%debug_str, logunit=diag_cs%chksum_iounit) - else if (is_stat) then + elseif (is_stat) then used = send_data(diag%fms_diag_id, locfield) elseif (diag_cs%ave_enabled) then used = send_data(diag%fms_diag_id, locfield, diag_cs%time_end) @@ -1333,7 +1336,7 @@ subroutine post_data_1d_k(diag_field_id, field, diag_cs, is_static) if (diag_cs%diag_as_chksum) then call zchksum(locfield, diag%debug_str, logunit=diag_cs%chksum_iounit) - else if (is_stat) then + elseif (is_stat) then used = send_data(diag%fms_diag_id, locfield) elseif (diag_cs%ave_enabled) then used = send_data(diag%fms_diag_id, locfield, diag_cs%time_end, weight=diag_cs%time_int) @@ -1448,38 +1451,38 @@ subroutine post_data_2d_low(diag, field, diag_cs, is_static, mask) endif if (present(mask)) then - locmask => mask - elseif(.NOT. is_stat) then - if(associated(diag%axes%mask2d)) locmask => diag%axes%mask2d + locmask => mask + elseif (.NOT. is_stat) then + if (associated(diag%axes%mask2d)) locmask => diag%axes%mask2d endif dl=1 - if(.NOT. is_stat) dl = diag%axes%downsample_level !static field downsample i not supported yet + if (.NOT. is_stat) dl = diag%axes%downsample_level !static field downsample i not supported yet !Downsample the diag field and mask (if present) if (dl > 1) then - isv_o=isv ; jsv_o=jsv - call downsample_diag_field(locfield, locfield_dsamp, dl, diag_cs, diag,isv,iev,jsv,jev, mask) - if ((diag%conversion_factor /= 0.) .and. (diag%conversion_factor /= 1.)) deallocate( locfield ) - locfield => locfield_dsamp - if (present(mask)) then - call downsample_field_2d(locmask, locmask_dsamp, dl, MSK, locmask, diag_cs,diag,isv_o,jsv_o,isv,iev,jsv,jev) - locmask => locmask_dsamp - elseif(associated(diag%axes%dsamp(dl)%mask2d)) then - locmask => diag%axes%dsamp(dl)%mask2d - endif + isv_o = isv ; jsv_o = jsv + call downsample_diag_field(locfield, locfield_dsamp, dl, diag_cs, diag,isv,iev,jsv,jev, mask) + if ((diag%conversion_factor /= 0.) .and. (diag%conversion_factor /= 1.)) deallocate( locfield ) + locfield => locfield_dsamp + if (present(mask)) then + call downsample_field_2d(locmask, locmask_dsamp, dl, MSK, locmask, diag_cs,diag,isv_o,jsv_o,isv,iev,jsv,jev) + locmask => locmask_dsamp + elseif (associated(diag%axes%dsamp(dl)%mask2d)) then + locmask => diag%axes%dsamp(dl)%mask2d + endif endif if (diag_cs%diag_as_chksum) then if (diag%axes%is_h_point) then call hchksum(locfield, diag%debug_str, diag_cs%G%HI, & logunit=diag_cs%chksum_iounit) - else if (diag%axes%is_u_point) then + elseif (diag%axes%is_u_point) then call uchksum(locfield, diag%debug_str, diag_cs%G%HI, & logunit=diag_cs%chksum_iounit) - else if (diag%axes%is_v_point) then + elseif (diag%axes%is_v_point) then call vchksum(locfield, diag%debug_str, diag_cs%G%HI, & logunit=diag_cs%chksum_iounit) - else if (diag%axes%is_q_point) then + elseif (diag%axes%is_q_point) then call Bchksum(locfield, diag%debug_str, diag_cs%G%HI, & logunit=diag_cs%chksum_iounit) else @@ -1735,25 +1738,25 @@ subroutine post_data_3d_low(diag, field, diag_cs, is_static, mask) endif if (present(mask)) then - locmask => mask - elseif(associated(diag%axes%mask3d)) then - locmask => diag%axes%mask3d + locmask => mask + elseif (associated(diag%axes%mask3d)) then + locmask => diag%axes%mask3d endif dl=1 - if(.NOT. is_stat) dl = diag%axes%downsample_level !static field downsample i not supported yet + if (.NOT. is_stat) dl = diag%axes%downsample_level !static field downsample i not supported yet !Downsample the diag field and mask (if present) if (dl > 1) then - isv_o=isv ; jsv_o=jsv - call downsample_diag_field(locfield, locfield_dsamp, dl, diag_cs, diag,isv,iev,jsv,jev, mask) - if ((diag%conversion_factor /= 0.) .and. (diag%conversion_factor /= 1.)) deallocate( locfield ) - locfield => locfield_dsamp - if (present(mask)) then - call downsample_field_3d(locmask, locmask_dsamp, dl, MSK, locmask, diag_cs,diag,isv_o,jsv_o,isv,iev,jsv,jev) - locmask => locmask_dsamp - elseif(associated(diag%axes%dsamp(dl)%mask3d)) then - locmask => diag%axes%dsamp(dl)%mask3d - endif + isv_o = isv ; jsv_o = jsv + call downsample_diag_field(locfield, locfield_dsamp, dl, diag_cs, diag,isv,iev,jsv,jev, mask) + if ((diag%conversion_factor /= 0.) .and. (diag%conversion_factor /= 1.)) deallocate( locfield ) + locfield => locfield_dsamp + if (present(mask)) then + call downsample_field_3d(locmask, locmask_dsamp, dl, MSK, locmask, diag_cs,diag,isv_o,jsv_o,isv,iev,jsv,jev) + locmask => locmask_dsamp + elseif (associated(diag%axes%dsamp(dl)%mask3d)) then + locmask => diag%axes%dsamp(dl)%mask3d + endif endif if (diag%fms_diag_id>0) then @@ -1761,13 +1764,13 @@ subroutine post_data_3d_low(diag, field, diag_cs, is_static, mask) if (diag%axes%is_h_point) then call hchksum(locfield, diag%debug_str, diag_cs%G%HI, & logunit=diag_cs%chksum_iounit) - else if (diag%axes%is_u_point) then + elseif (diag%axes%is_u_point) then call uchksum(locfield, diag%debug_str, diag_cs%G%HI, & logunit=diag_cs%chksum_iounit) - else if (diag%axes%is_v_point) then + elseif (diag%axes%is_v_point) then call vchksum(locfield, diag%debug_str, diag_cs%G%HI, & logunit=diag_cs%chksum_iounit) - else if (diag%axes%is_q_point) then + elseif (diag%axes%is_q_point) then call Bchksum(locfield, diag%debug_str, diag_cs%G%HI, & logunit=diag_cs%chksum_iounit) else @@ -1941,12 +1944,12 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time verbose, do_not_log, err_msg, interp_method, tile_count, cmor_field_name, & cmor_long_name, cmor_units, cmor_standard_name, cell_methods, & x_cell_method, y_cell_method, v_cell_method, conversion, v_extensive) - character(len=*), intent(in) :: module_name !< Name of this module, usually "ocean_model" - !! or "ice_shelf_model" - character(len=*), intent(in) :: field_name !< Name of the diagnostic field - type(axes_grp), target, intent(in) :: axes_in !< Container w/ up to 3 integer handles that - !! indicates axes for this field - type(time_type), intent(in) :: init_time !< Time at which a field is first available? + character(len=*), intent(in) :: module_name !< Name of this module, usually "ocean_model" + !! or "ice_shelf_model" + character(len=*), intent(in) :: field_name !< Name of the diagnostic field + type(axes_grp), target, intent(in) :: axes_in !< Container w/ up to 3 integer handles that + !! indicates axes for this field + type(time_type), intent(in) :: init_time !< Time at which a field is first available? character(len=*), optional, intent(in) :: long_name !< Long name of a field. character(len=*), optional, intent(in) :: units !< Units of a field. character(len=*), optional, intent(in) :: standard_name !< Standardized name associated with a field @@ -1984,7 +1987,10 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time type(axes_grp), pointer :: remap_axes => null() type(axes_grp), pointer :: axes => null() integer :: dm_id, i, dl + character(len=256) :: msg, cm_string character(len=256) :: new_module_name + character(len=480) :: module_list, var_list, variants + integer :: num_modnm, num_varnm logical :: active axes => axes_in @@ -1995,23 +2001,26 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time dm_id = -1 if (axes_in%id == diag_cs%axesTL%id) then - axes => diag_cs%axesTL + axes => diag_cs%axesTL elseif (axes_in%id == diag_cs%axesBL%id) then - axes => diag_cs%axesBL + axes => diag_cs%axesBL elseif (axes_in%id == diag_cs%axesCuL%id ) then - axes => diag_cs%axesCuL + axes => diag_cs%axesCuL elseif (axes_in%id == diag_cs%axesCvL%id) then - axes => diag_cs%axesCvL + axes => diag_cs%axesCvL elseif (axes_in%id == diag_cs%axesTi%id) then - axes => diag_cs%axesTi + axes => diag_cs%axesTi elseif (axes_in%id == diag_cs%axesBi%id) then - axes => diag_cs%axesBi + axes => diag_cs%axesBi elseif (axes_in%id == diag_cs%axesCui%id ) then - axes => diag_cs%axesCui + axes => diag_cs%axesCui elseif (axes_in%id == diag_cs%axesCvi%id) then - axes => diag_cs%axesCvi + axes => diag_cs%axesCvi endif + module_list = "{$M" + num_modnm = 1 + ! Register the native diagnostic active = register_diag_field_expand_cmor(dm_id, module_name, field_name, axes, & init_time, long_name=long_name, units=units, missing_value=MOM_missing_value, & @@ -2023,6 +2032,21 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time cell_methods=cell_methods, x_cell_method=x_cell_method, & y_cell_method=y_cell_method, v_cell_method=v_cell_method, & conversion=conversion, v_extensive=v_extensive) + if (associated(axes%xyave_axes)) then + num_varnm = 2 ; var_list = "{$V,$V_xyave" + else + num_varnm = 1 ; var_list = "{$V" + endif + if (present(cmor_field_name)) then + if (associated(axes%xyave_axes)) then + num_varnm = num_varnm + 2 + var_list = trim(var_list)//","//trim(cmor_field_name)//","//trim(cmor_field_name)//"_xyave" + else + num_varnm = num_varnm + 1 + var_list = trim(var_list)//","//trim(cmor_field_name) + endif + endif + var_list = trim(var_list)//"}" ! For each diagnostic coordinate register the diagnostic again under a different module name do i=1,diag_cs%num_diag_coords @@ -2032,21 +2056,21 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time if (axes_in%rank == 3) then remap_axes => null() if ((axes_in%id == diag_cs%axesTL%id)) then - remap_axes => diag_cs%remap_axesTL(i) + remap_axes => diag_cs%remap_axesTL(i) elseif (axes_in%id == diag_cs%axesBL%id) then - remap_axes => diag_cs%remap_axesBL(i) + remap_axes => diag_cs%remap_axesBL(i) elseif (axes_in%id == diag_cs%axesCuL%id ) then - remap_axes => diag_cs%remap_axesCuL(i) + remap_axes => diag_cs%remap_axesCuL(i) elseif (axes_in%id == diag_cs%axesCvL%id) then - remap_axes => diag_cs%remap_axesCvL(i) + remap_axes => diag_cs%remap_axesCvL(i) elseif (axes_in%id == diag_cs%axesTi%id) then - remap_axes => diag_cs%remap_axesTi(i) + remap_axes => diag_cs%remap_axesTi(i) elseif (axes_in%id == diag_cs%axesBi%id) then - remap_axes => diag_cs%remap_axesBi(i) + remap_axes => diag_cs%remap_axesBi(i) elseif (axes_in%id == diag_cs%axesCui%id ) then - remap_axes => diag_cs%remap_axesCui(i) + remap_axes => diag_cs%remap_axesCui(i) elseif (axes_in%id == diag_cs%axesCvi%id) then - remap_axes => diag_cs%remap_axesCvi(i) + remap_axes => diag_cs%remap_axesCvi(i) endif ! When the MOM_diag_to_Z module has been obsoleted we can assume remap_axes will ! always exist but in the mean-time we have to do this check: @@ -2066,6 +2090,9 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time if (active) then call diag_remap_set_active(diag_cs%diag_remap_cs(i)) endif + ! module_list = trim(module_list)//","//trim(new_module_name) + module_list = trim(module_list)//",$M_"//trim(diag_cs%diag_remap_cs(i)%diag_module_suffix) + num_modnm = num_modnm + 1 endif ! remap_axes%needs_remapping endif ! associated(remap_axes) endif ! axes%rank == 3 @@ -2073,46 +2100,46 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time !Register downsampled diagnostics do dl=2,MAX_DSAMP_LEV - ! Do not attempt to checksum the downsampled diagnostics - if (diag_cs%diag_as_chksum) cycle + ! Do not attempt to checksum the downsampled diagnostics + if (diag_cs%diag_as_chksum) cycle - new_module_name = trim(module_name)//'_d2' + new_module_name = trim(module_name)//'_d2' - if (axes_in%rank == 3 .or. axes_in%rank == 2 ) then - axes => null() - if (axes_in%id == diag_cs%axesTL%id) then - axes => diag_cs%dsamp(dl)%axesTL - elseif (axes_in%id == diag_cs%axesBL%id) then - axes => diag_cs%dsamp(dl)%axesBL - elseif (axes_in%id == diag_cs%axesCuL%id ) then - axes => diag_cs%dsamp(dl)%axesCuL - elseif (axes_in%id == diag_cs%axesCvL%id) then - axes => diag_cs%dsamp(dl)%axesCvL - elseif (axes_in%id == diag_cs%axesTi%id) then - axes => diag_cs%dsamp(dl)%axesTi - elseif (axes_in%id == diag_cs%axesBi%id) then - axes => diag_cs%dsamp(dl)%axesBi - elseif (axes_in%id == diag_cs%axesCui%id ) then - axes => diag_cs%dsamp(dl)%axesCui - elseif (axes_in%id == diag_cs%axesCvi%id) then - axes => diag_cs%dsamp(dl)%axesCvi - elseif (axes_in%id == diag_cs%axesT1%id) then - axes => diag_cs%dsamp(dl)%axesT1 - elseif (axes_in%id == diag_cs%axesB1%id) then - axes => diag_cs%dsamp(dl)%axesB1 - elseif (axes_in%id == diag_cs%axesCu1%id ) then - axes => diag_cs%dsamp(dl)%axesCu1 - elseif (axes_in%id == diag_cs%axesCv1%id) then - axes => diag_cs%dsamp(dl)%axesCv1 - else - !Niki: Should we worry about these, e.g., diag_to_Z_CS? - call MOM_error(WARNING,"register_diag_field: Could not find a proper axes for " & - //trim( new_module_name)//"-"//trim(field_name)) - endif - endif - ! Register the native diagnostic - if (associated(axes)) then - active = register_diag_field_expand_cmor(dm_id, new_module_name, field_name, axes, & + if (axes_in%rank == 3 .or. axes_in%rank == 2 ) then + axes => null() + if (axes_in%id == diag_cs%axesTL%id) then + axes => diag_cs%dsamp(dl)%axesTL + elseif (axes_in%id == diag_cs%axesBL%id) then + axes => diag_cs%dsamp(dl)%axesBL + elseif (axes_in%id == diag_cs%axesCuL%id ) then + axes => diag_cs%dsamp(dl)%axesCuL + elseif (axes_in%id == diag_cs%axesCvL%id) then + axes => diag_cs%dsamp(dl)%axesCvL + elseif (axes_in%id == diag_cs%axesTi%id) then + axes => diag_cs%dsamp(dl)%axesTi + elseif (axes_in%id == diag_cs%axesBi%id) then + axes => diag_cs%dsamp(dl)%axesBi + elseif (axes_in%id == diag_cs%axesCui%id ) then + axes => diag_cs%dsamp(dl)%axesCui + elseif (axes_in%id == diag_cs%axesCvi%id) then + axes => diag_cs%dsamp(dl)%axesCvi + elseif (axes_in%id == diag_cs%axesT1%id) then + axes => diag_cs%dsamp(dl)%axesT1 + elseif (axes_in%id == diag_cs%axesB1%id) then + axes => diag_cs%dsamp(dl)%axesB1 + elseif (axes_in%id == diag_cs%axesCu1%id ) then + axes => diag_cs%dsamp(dl)%axesCu1 + elseif (axes_in%id == diag_cs%axesCv1%id) then + axes => diag_cs%dsamp(dl)%axesCv1 + else + !Niki: Should we worry about these, e.g., diag_to_Z_CS? + call MOM_error(WARNING,"register_diag_field: Could not find a proper axes for " & + //trim(new_module_name)//"-"//trim(field_name)) + endif + endif + ! Register the native diagnostic + if (associated(axes)) then + active = register_diag_field_expand_cmor(dm_id, new_module_name, field_name, axes, & init_time, long_name=long_name, units=units, missing_value=MOM_missing_value, & range=range, mask_variant=mask_variant, standard_name=standard_name, & verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, & @@ -2122,57 +2149,83 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time cell_methods=cell_methods, x_cell_method=x_cell_method, & y_cell_method=y_cell_method, v_cell_method=v_cell_method, & conversion=conversion, v_extensive=v_extensive) - endif + ! module_list = trim(module_list)//","//trim(new_module_name) + module_list = trim(module_list)//",$M_2d" + num_modnm = num_modnm + 1 + endif - ! For each diagnostic coordinate register the diagnostic again under a different module name - do i=1,diag_cs%num_diag_coords - new_module_name = trim(module_name)//'_'//trim(diag_cs%diag_remap_cs(i)%diag_module_suffix)//'_d2' - - ! Register diagnostics remapped to z vertical coordinate - if (axes_in%rank == 3) then - remap_axes => null() - if ((axes_in%id == diag_cs%axesTL%id)) then - remap_axes => diag_cs%dsamp(dl)%remap_axesTL(i) - elseif (axes_in%id == diag_cs%axesBL%id) then - remap_axes => diag_cs%dsamp(dl)%remap_axesBL(i) - elseif (axes_in%id == diag_cs%axesCuL%id ) then - remap_axes => diag_cs%dsamp(dl)%remap_axesCuL(i) - elseif (axes_in%id == diag_cs%axesCvL%id) then - remap_axes => diag_cs%dsamp(dl)%remap_axesCvL(i) - elseif (axes_in%id == diag_cs%axesTi%id) then - remap_axes => diag_cs%dsamp(dl)%remap_axesTi(i) - elseif (axes_in%id == diag_cs%axesBi%id) then - remap_axes => diag_cs%dsamp(dl)%remap_axesBi(i) - elseif (axes_in%id == diag_cs%axesCui%id ) then - remap_axes => diag_cs%dsamp(dl)%remap_axesCui(i) - elseif (axes_in%id == diag_cs%axesCvi%id) then - remap_axes => diag_cs%dsamp(dl)%remap_axesCvi(i) - endif - - ! When the MOM_diag_to_Z module has been obsoleted we can assume remap_axes will - ! always exist but in the mean-time we have to do this check: - ! call assert(associated(remap_axes), 'register_diag_field: remap_axes not set') - if (associated(remap_axes)) then - if (remap_axes%needs_remapping .or. remap_axes%needs_interpolating) then - active = register_diag_field_expand_cmor(dm_id, new_module_name, field_name, remap_axes, & - init_time, long_name=long_name, units=units, missing_value=MOM_missing_value, & - range=range, mask_variant=mask_variant, standard_name=standard_name, & - verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, & - interp_method=interp_method, tile_count=tile_count, & - cmor_field_name=cmor_field_name, cmor_long_name=cmor_long_name, & - cmor_units=cmor_units, cmor_standard_name=cmor_standard_name, & - cell_methods=cell_methods, x_cell_method=x_cell_method, & - y_cell_method=y_cell_method, v_cell_method=v_cell_method, & - conversion=conversion, v_extensive=v_extensive) - if (active) then - call diag_remap_set_active(diag_cs%diag_remap_cs(i)) - endif - endif ! remap_axes%needs_remapping - endif ! associated(remap_axes) - endif ! axes%rank == 3 - enddo ! i + ! For each diagnostic coordinate register the diagnostic again under a different module name + do i=1,diag_cs%num_diag_coords + new_module_name = trim(module_name)//'_'//trim(diag_cs%diag_remap_cs(i)%diag_module_suffix)//'_d2' + + ! Register diagnostics remapped to z vertical coordinate + if (axes_in%rank == 3) then + remap_axes => null() + if ((axes_in%id == diag_cs%axesTL%id)) then + remap_axes => diag_cs%dsamp(dl)%remap_axesTL(i) + elseif (axes_in%id == diag_cs%axesBL%id) then + remap_axes => diag_cs%dsamp(dl)%remap_axesBL(i) + elseif (axes_in%id == diag_cs%axesCuL%id ) then + remap_axes => diag_cs%dsamp(dl)%remap_axesCuL(i) + elseif (axes_in%id == diag_cs%axesCvL%id) then + remap_axes => diag_cs%dsamp(dl)%remap_axesCvL(i) + elseif (axes_in%id == diag_cs%axesTi%id) then + remap_axes => diag_cs%dsamp(dl)%remap_axesTi(i) + elseif (axes_in%id == diag_cs%axesBi%id) then + remap_axes => diag_cs%dsamp(dl)%remap_axesBi(i) + elseif (axes_in%id == diag_cs%axesCui%id ) then + remap_axes => diag_cs%dsamp(dl)%remap_axesCui(i) + elseif (axes_in%id == diag_cs%axesCvi%id) then + remap_axes => diag_cs%dsamp(dl)%remap_axesCvi(i) + endif + + ! When the MOM_diag_to_Z module has been obsoleted we can assume remap_axes will + ! always exist but in the mean-time we have to do this check: + ! call assert(associated(remap_axes), 'register_diag_field: remap_axes not set') + if (associated(remap_axes)) then + if (remap_axes%needs_remapping .or. remap_axes%needs_interpolating) then + active = register_diag_field_expand_cmor(dm_id, new_module_name, field_name, remap_axes, & + init_time, long_name=long_name, units=units, missing_value=MOM_missing_value, & + range=range, mask_variant=mask_variant, standard_name=standard_name, & + verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, & + interp_method=interp_method, tile_count=tile_count, & + cmor_field_name=cmor_field_name, cmor_long_name=cmor_long_name, & + cmor_units=cmor_units, cmor_standard_name=cmor_standard_name, & + cell_methods=cell_methods, x_cell_method=x_cell_method, & + y_cell_method=y_cell_method, v_cell_method=v_cell_method, & + conversion=conversion, v_extensive=v_extensive) + if (active) then + call diag_remap_set_active(diag_cs%diag_remap_cs(i)) + endif + ! module_list = trim(module_list)//","//trim(new_module_name) + module_list = trim(module_list)//",$M_"//trim(diag_cs%diag_remap_cs(i)%diag_module_suffix)//"_d2" + num_modnm = num_modnm + 1 + endif ! remap_axes%needs_remapping + endif ! associated(remap_axes) + endif ! axes%rank == 3 + enddo ! i enddo + if (is_root_pe() .and. (diag_CS%available_diag_doc_unit > 0) .and. .not.diag_CS%full_diag_list) then + msg = '' + if (present(cmor_field_name)) msg = 'CMOR equivalent is "'//trim(cmor_field_name)//'"' + call attach_cell_methods(-1, axes, cm_string, cell_methods, & + x_cell_method, y_cell_method, v_cell_method, & + v_extensive=v_extensive) + module_list = trim(module_list)//"}" + if ((num_modnm > 1) .and. (num_varnm > 1)) then + variants = "mod:"//trim(module_list)//", var:"//trim(var_list) + elseif (num_modnm > 1) then + variants = "mod:"//trim(module_list) + elseif (num_varnm > 1) then + variants = "var:"//trim(var_list) + else + variants = "" + endif + call log_available_diag(dm_id>0, module_name, field_name, cm_string, msg, diag_CS, & + long_name, units, standard_name, variants=variants) + endif + register_diag_field = dm_id end function register_diag_field @@ -2244,7 +2297,7 @@ logical function register_diag_field_expand_cmor(dm_id, module_name, field_name, call attach_cell_methods(fms_id, axes, cm_string, cell_methods, & x_cell_method, y_cell_method, v_cell_method, & v_extensive=v_extensive) - if (is_root_pe() .and. diag_CS%available_diag_doc_unit > 0) then + if (is_root_pe() .and. (diag_CS%available_diag_doc_unit > 0) .and. diag_CS%full_diag_list) then msg = '' if (present(cmor_field_name)) msg = 'CMOR equivalent is "'//trim(cmor_field_name)//'"' call log_available_diag(fms_id>0, module_name, field_name, cm_string, & @@ -2262,7 +2315,7 @@ logical function register_diag_field_expand_cmor(dm_id, module_name, field_name, if (.not. diag_cs%diag_as_chksum) & call attach_cell_methods(fms_xyave_id, axes%xyave_axes, cm_string, & cell_methods, v_cell_method, v_extensive=v_extensive) - if (is_root_pe() .and. diag_CS%available_diag_doc_unit > 0) then + if (is_root_pe() .and. (diag_CS%available_diag_doc_unit > 0) .and. diag_CS%full_diag_list) then msg = '' if (present(cmor_field_name)) msg = 'CMOR equivalent is "'//trim(cmor_field_name)//'_xyave"' call log_available_diag(fms_xyave_id>0, module_name, trim(field_name)//'_xyave', cm_string, & @@ -2306,7 +2359,7 @@ logical function register_diag_field_expand_cmor(dm_id, module_name, field_name, call attach_cell_methods(fms_id, axes, cm_string, & cell_methods, x_cell_method, y_cell_method, v_cell_method, & v_extensive=v_extensive) - if (is_root_pe() .and. diag_CS%available_diag_doc_unit > 0) then + if (is_root_pe() .and. (diag_CS%available_diag_doc_unit > 0) .and. diag_CS%full_diag_list) then msg = 'native name is "'//trim(field_name)//'"' call log_available_diag(fms_id>0, module_name, cmor_field_name, cm_string, & msg, diag_CS, posted_cmor_long_name, posted_cmor_units, & @@ -2323,7 +2376,7 @@ logical function register_diag_field_expand_cmor(dm_id, module_name, field_name, err_msg=err_msg, interp_method=interp_method, tile_count=tile_count) call attach_cell_methods(fms_xyave_id, axes%xyave_axes, cm_string, & cell_methods, v_cell_method, v_extensive=v_extensive) - if (is_root_pe() .and. diag_CS%available_diag_doc_unit > 0) then + if (is_root_pe() .and. (diag_CS%available_diag_doc_unit > 0) .and. diag_CS%full_diag_list) then msg = 'native name is "'//trim(field_name)//'_xyave"' call log_available_diag(fms_xyave_id>0, module_name, trim(cmor_field_name)//'_xyave', & cm_string, msg, diag_CS, posted_cmor_long_name, posted_cmor_units, & @@ -2381,7 +2434,7 @@ integer function register_diag_field_expand_axes(module_name, field_name, axes, if (axes%diag_cs%diag_as_chksum) then fms_id = axes%diag_cs%num_chksum_diags + 1 axes%diag_cs%num_chksum_diags = fms_id - else if (present(interp_method) .or. axes%is_h_point) then + elseif (present(interp_method) .or. axes%is_h_point) then ! If interp_method is provided we must use it if (area_id>0) then if (volume_id>0) then @@ -2504,7 +2557,7 @@ subroutine add_xyz_method(diag, axes, x_cell_method, y_cell_method, v_cell_metho if (present(v_extensive)) then if (present(v_cell_method)) call MOM_error(FATAL, "attach_cell_methods: " // & 'Vertical cell method was specified along with the vertically extensive flag.') - if(v_extensive) then + if (v_extensive) then mstr='sum' else mstr='mean' @@ -2624,7 +2677,7 @@ subroutine attach_cell_methods(id, axes, ostring, cell_methods, & ostring = trim(adjustl(ostring))//' '//trim(axis_name)//':'//trim(v_cell_method) endif elseif (present(v_extensive)) then - if(v_extensive) then + if (v_extensive) then if (axes%rank==1) then call get_diag_axis_name(axes%handles(1), axis_name) elseif (axes%rank==3) then @@ -2744,9 +2797,15 @@ function register_scalar_field(module_name, field_name, init_time, diag_cs, & ! Document diagnostics in list of available diagnostics if (is_root_pe() .and. diag_CS%available_diag_doc_unit > 0) then - call log_available_diag(associated(diag), module_name, field_name, '', '', diag_CS, & - long_name, units, standard_name) - if (present(cmor_field_name)) then + if (present(cmor_field_name) .and. .not.diag_CS%full_diag_list) then + call log_available_diag(associated(diag), module_name, field_name, '', '', diag_CS, & + long_name, units, standard_name, & + variants="var:{$V,"//trim(cmor_field_name)//"}") + else + call log_available_diag(associated(diag), module_name, field_name, '', '', diag_CS, & + long_name, units, standard_name) + endif + if (present(cmor_field_name) .and. diag_CS%full_diag_list) then call log_available_diag(associated(cmor_diag), module_name, cmor_field_name, & '', '', diag_CS, posted_cmor_long_name, posted_cmor_units, & posted_cmor_standard_name) @@ -2891,9 +2950,15 @@ function register_static_field(module_name, field_name, axes, & ! Document diagnostics in list of available diagnostics if (is_root_pe() .and. diag_CS%available_diag_doc_unit > 0) then - call log_available_diag(associated(diag), module_name, field_name, '', '', diag_CS, & - long_name, units, standard_name) - if (present(cmor_field_name)) then + if (present(cmor_field_name) .and. .not.diag_CS%full_diag_list) then + call log_available_diag(associated(diag), module_name, field_name, '', '', diag_CS, & + long_name, units, standard_name, & + variants="var:{$V,"//trim(cmor_field_name)//"}") + else + call log_available_diag(associated(diag), module_name, field_name, '', '', diag_CS, & + long_name, units, standard_name) + endif + if (present(cmor_field_name) .and. diag_CS%full_diag_list) then call log_available_diag(associated(cmor_diag), module_name, cmor_field_name, & '', '', diag_CS, posted_cmor_long_name, posted_cmor_units, & posted_cmor_standard_name) @@ -2910,7 +2975,7 @@ subroutine describe_option(opt_name, value, diag_CS) character(len=*), intent(in) :: value !< A character string with the setting of the option. type(diag_ctrl), intent(in) :: diag_CS !< Structure used to regulate diagnostic output - character(len=240) :: mesg + character(len=480) :: mesg integer :: len_ind len_ind = len_trim(value) ! Add error handling for long values? @@ -2944,75 +3009,46 @@ function ocean_register_diag(var_desc, G, diag_CS, day) case ("L") select case (hor_grid) - case ("q") - axes => diag_cs%axesBL - case ("h") - axes => diag_cs%axesTL - case ("u") - axes => diag_cs%axesCuL - case ("v") - axes => diag_cs%axesCvL - case ("Bu") - axes => diag_cs%axesBL - case ("T") - axes => diag_cs%axesTL - case ("Cu") - axes => diag_cs%axesCuL - case ("Cv") - axes => diag_cs%axesCvL - case ("z") - axes => diag_cs%axeszL - case default - call MOM_error(FATAL, "ocean_register_diag: " // & - "unknown hor_grid component "//trim(hor_grid)) + case ("q") ; axes => diag_cs%axesBL + case ("h") ; axes => diag_cs%axesTL + case ("u") ; axes => diag_cs%axesCuL + case ("v") ; axes => diag_cs%axesCvL + case ("Bu") ; axes => diag_cs%axesBL + case ("T") ; axes => diag_cs%axesTL + case ("Cu") ; axes => diag_cs%axesCuL + case ("Cv") ; axes => diag_cs%axesCvL + case ("z") ; axes => diag_cs%axeszL + case default ; call MOM_error(FATAL, "ocean_register_diag: " // & + "unknown hor_grid component "//trim(hor_grid)) end select case ("i") select case (hor_grid) - case ("q") - axes => diag_cs%axesBi - case ("h") - axes => diag_cs%axesTi - case ("u") - axes => diag_cs%axesCui - case ("v") - axes => diag_cs%axesCvi - case ("Bu") - axes => diag_cs%axesBi - case ("T") - axes => diag_cs%axesTi - case ("Cu") - axes => diag_cs%axesCui - case ("Cv") - axes => diag_cs%axesCvi - case ("z") - axes => diag_cs%axeszi - case default - call MOM_error(FATAL, "ocean_register_diag: " // & - "unknown hor_grid component "//trim(hor_grid)) + case ("q") ; axes => diag_cs%axesBi + case ("h") ; axes => diag_cs%axesTi + case ("u") ; axes => diag_cs%axesCui + case ("v") ; axes => diag_cs%axesCvi + case ("Bu") ; axes => diag_cs%axesBi + case ("T") ; axes => diag_cs%axesTi + case ("Cu") ; axes => diag_cs%axesCui + case ("Cv") ; axes => diag_cs%axesCvi + case ("z") ; axes => diag_cs%axeszi + case default ; call MOM_error(FATAL, "ocean_register_diag: " // & + "unknown hor_grid component "//trim(hor_grid)) end select case ("1") select case (hor_grid) - case ("q") - axes => diag_cs%axesB1 - case ("h") - axes => diag_cs%axesT1 - case ("u") - axes => diag_cs%axesCu1 - case ("v") - axes => diag_cs%axesCv1 - case ("Bu") - axes => diag_cs%axesB1 - case ("T") - axes => diag_cs%axesT1 - case ("Cu") - axes => diag_cs%axesCu1 - case ("Cv") - axes => diag_cs%axesCv1 - case default - call MOM_error(FATAL, "ocean_register_diag: " // & - "unknown hor_grid component "//trim(hor_grid)) + case ("q") ; axes => diag_cs%axesB1 + case ("h") ; axes => diag_cs%axesT1 + case ("u") ; axes => diag_cs%axesCu1 + case ("v") ; axes => diag_cs%axesCv1 + case ("Bu") ; axes => diag_cs%axesB1 + case ("T") ; axes => diag_cs%axesT1 + case ("Cu") ; axes => diag_cs%axesCu1 + case ("Cv") ; axes => diag_cs%axesCv1 + case default ; call MOM_error(FATAL, "ocean_register_diag: " // & + "unknown hor_grid component "//trim(hor_grid)) end select case default @@ -3190,6 +3226,12 @@ subroutine diag_mediator_init(G, GV, US, nz, param_file, diag_cs, doc_file_dir) if ((.not.opened) .or. (ios /= 0)) then call MOM_error(FATAL, "Failed to open available diags file "//trim(doc_path)//".") endif + + call get_param(param_file, mdl, "FULL_AVAILABLE_DIAGS_LIST", diag_CS%full_diag_list, & + "If true, write entries for every variant of a diagnostic to the available //"& + "diags file; otherwise only write an entry for the primary diagnostics."//& + "ocean diagnostics that can be included in a diag_table.", & + default=.true.) endif endif @@ -3540,7 +3582,7 @@ end subroutine alloc_diag_with_id !> Log a diagnostic to the available diagnostics file. subroutine log_available_diag(used, module_name, field_name, cell_methods_string, comment, & - diag_CS, long_name, units, standard_name) + diag_CS, long_name, units, standard_name, variants) logical, intent(in) :: used !< Whether this diagnostic was in the diag_table or not character(len=*), intent(in) :: module_name !< Name of the diagnostic module character(len=*), intent(in) :: field_name !< Name of this diagnostic field @@ -3550,6 +3592,8 @@ subroutine log_available_diag(used, module_name, field_name, cell_methods_string character(len=*), optional, intent(in) :: long_name !< CF long name of diagnostic character(len=*), optional, intent(in) :: units !< Units for diagnostic character(len=*), optional, intent(in) :: standard_name !< CF standardized name of diagnostic + character(len=*), optional, intent(in) :: variants !< Alternate modules and variable names for + !! this diagnostic and derived diagnostics ! Local variables character(len=240) :: mesg @@ -3569,7 +3613,9 @@ subroutine log_available_diag(used, module_name, field_name, cell_methods_string call describe_option("standard_name", standard_name, diag_CS) if (len(trim((cell_methods_string)))>0) & call describe_option("cell_methods", trim(cell_methods_string), diag_CS) - + if (present(variants)) then ; if (len(trim(variants)) > 0) then + call describe_option("variants", variants, diag_CS) + endif ; endif end subroutine log_available_diag !> Log the diagnostic chksum to the chksum diag file @@ -3778,8 +3824,8 @@ subroutine downsample_diag_indices_get(fo1, fo2, dl, diag_cs, isv, iev, jsv, jev !We assume that the compute domain can be subdivided to dl*dl cells, hence avoiding the need of halo updates. !We want this check to error out only if there was a downsampled diagnostics requested and about to post that is !why the check is here and not in the init routines. This check need to be done only once, hence the outer if. - if(first_check) then - if(mod(diag_cs%ie-diag_cs%is+1, dl) .ne. 0 .OR. mod(diag_cs%je-diag_cs%js+1, dl) .ne. 0) then + if (first_check) then + if (mod(diag_cs%ie-diag_cs%is+1, dl) /= 0 .OR. mod(diag_cs%je-diag_cs%js+1, dl) /= 0) then write (mesg,*) "Non-commensurate downsampled domain is not supported. "//& "Please choose a layout such that NIGLOBAL/Layout_X and NJGLOBAL/Layout_Y are both divisible by dl=",dl,& " Current domain extents: ", diag_cs%is,diag_cs%ie, diag_cs%js,diag_cs%je @@ -3849,10 +3895,10 @@ subroutine downsample_diag_field_3d(locfield, locfield_dsamp, dl, diag_cs, diag, locmask => NULL() !Get the correct indices corresponding to input field !Shape of the input diag field - f1=size(locfield,1) - f2=size(locfield,2) + f1 = size(locfield,1) + f2 = size(locfield,2) !Save the extents of the original (fine) domain - isv_o=isv;jsv_o=jsv + isv_o = isv ; jsv_o = jsv !Get the shape of the downsampled field and overwrite isv,iev,jsv,jev with them call downsample_diag_indices_get(f1,f2, dl, diag_cs,isv,iev,jsv,jev) !Set the non-downsampled mask, it must be associated and initialized @@ -3890,19 +3936,19 @@ subroutine downsample_diag_field_2d(locfield, locfield_dsamp, dl, diag_cs, diag, locmask => NULL() !Get the correct indices corresponding to input field !Shape of the input diag field - f1=size(locfield,1) - f2=size(locfield,2) + f1 = size(locfield,1) + f2 = size(locfield,2) !Save the extents of the original (fine) domain - isv_o=isv;jsv_o=jsv + isv_o = isv ; jsv_o = jsv !Get the shape of the downsampled field and overwrite isv,iev,jsv,jev with them call downsample_diag_indices_get(f1,f2, dl, diag_cs,isv,iev,jsv,jev) !Set the non-downsampled mask, it must be associated and initialized if (present(mask)) then - locmask => mask + locmask => mask elseif (associated(diag%axes%mask2d)) then - locmask => diag%axes%mask2d + locmask => diag%axes%mask2d else - call MOM_error(FATAL, "downsample_diag_field_2d: Cannot downsample without a mask!!! ") + call MOM_error(FATAL, "downsample_diag_field_2d: Cannot downsample without a mask!!! ") endif call downsample_field(locfield, locfield_dsamp, dl, diag%xyz_method, locmask, diag_cs,diag, & @@ -3991,7 +4037,7 @@ subroutine downsample_field_3d(field_in, field_out, dl, method, mask, diag_cs, d ! Fill the down sampled field on the down sampled diagnostics (almost always compuate) domain if (method == MMM) then - do k= ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d + do k=ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d i0 = isv_o+dl*(i-isv_d) j0 = jsv_o+dl*(j-jsv_d) ave = 0.0 @@ -4001,22 +4047,22 @@ subroutine downsample_field_3d(field_in, field_out, dl, method, mask, diag_cs, d weight = mask(ii,jj,k) * diag_cs%G%areaT(ii,jj) * diag_cs%h(ii,jj,k) total_weight = total_weight + weight ave = ave+field_in(ii,jj,k) * weight - enddo; enddo + enddo ; enddo field_out(i,j,k) = ave/(total_weight + eps_vol) !Avoid zero mask at all aggregating cells where ave=0.0 - enddo; enddo; enddo + enddo ; enddo ; enddo elseif (method == SSS) then !e.g., volcello - do k= ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d + do k=ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d i0 = isv_o+dl*(i-isv_d) j0 = jsv_o+dl*(j-jsv_d) ave = 0.0 do jj=j0,j0+dl-1 ; do ii=i0,i0+dl-1 weight = mask(ii,jj,k) ave = ave+field_in(ii,jj,k)*weight - enddo; enddo + enddo ; enddo field_out(i,j,k) = ave !Masked Sum (total_weight=1) - enddo; enddo; enddo - elseif(method == MMP .or. method == MMS) then !e.g., T_advection_xy - do k= ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d + enddo ; enddo ; enddo + elseif (method == MMP .or. method == MMS) then !e.g., T_advection_xy + do k=ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d i0 = isv_o+dl*(i-isv_d) j0 = jsv_o+dl*(j-jsv_d) ave = 0.0 @@ -4026,49 +4072,49 @@ subroutine downsample_field_3d(field_in, field_out, dl, method, mask, diag_cs, d weight = mask(ii,jj,k) * diag_cs%G%areaT(ii,jj) total_weight = total_weight + weight ave = ave+field_in(ii,jj,k)*weight - enddo; enddo + enddo ; enddo field_out(i,j,k) = ave / (total_weight+eps_area) !Avoid zero mask at all aggregating cells where ave=0.0 - enddo; enddo; enddo - elseif(method == PMM) then - do k= ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d + enddo ; enddo ; enddo + elseif (method == PMM) then + do k=ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d i0 = isv_o+dl*(i-isv_d) j0 = jsv_o+dl*(j-jsv_d) ave = 0.0 total_weight = 0.0 ii=i0 do jj=j0,j0+dl-1 - weight =mask(ii,jj,k) * diag_cs%G%dyCu(ii,jj) * diag_cs%h(ii,jj,k) + weight = mask(ii,jj,k) * diag_cs%G%dyCu(ii,jj) * diag_cs%h(ii,jj,k) total_weight = total_weight +weight - ave=ave+field_in(ii,jj,k)*weight + ave = ave+field_in(ii,jj,k)*weight enddo field_out(i,j,k) = ave/(total_weight+eps_face) !Avoid zero mask at all aggregating cells where ave=0.0 - enddo; enddo; enddo - elseif(method == PSS) then !e.g. umo - do k= ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d + enddo ; enddo ; enddo + elseif (method == PSS) then !e.g. umo + do k=ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d i0 = isv_o+dl*(i-isv_d) j0 = jsv_o+dl*(j-jsv_d) ave = 0.0 ii=i0 do jj=j0,j0+dl-1 - weight =mask(ii,jj,k) - ave=ave+field_in(ii,jj,k)*weight + weight = mask(ii,jj,k) + ave = ave+field_in(ii,jj,k)*weight enddo field_out(i,j,k) = ave !Masked Sum (total_weight=1) - enddo; enddo; enddo - elseif(method == SPS) then !e.g. vmo - do k= ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d + enddo ; enddo ; enddo + elseif (method == SPS) then !e.g. vmo + do k=ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d i0 = isv_o+dl*(i-isv_d) j0 = jsv_o+dl*(j-jsv_d) ave = 0.0 jj=j0 do ii=i0,i0+dl-1 - weight =mask(ii,jj,k) - ave=ave+field_in(ii,jj,k)*weight + weight = mask(ii,jj,k) + ave = ave+field_in(ii,jj,k)*weight enddo field_out(i,j,k) = ave !Masked Sum (total_weight=1) - enddo; enddo; enddo - elseif(method == MPM) then - do k= ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d + enddo ; enddo ; enddo + elseif (method == MPM) then + do k=ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d i0 = isv_o+dl*(i-isv_d) j0 = jsv_o+dl*(j-jsv_d) ave = 0.0 @@ -4077,21 +4123,21 @@ subroutine downsample_field_3d(field_in, field_out, dl, method, mask, diag_cs, d do ii=i0,i0+dl-1 weight = mask(ii,jj,k) * diag_cs%G%dxCv(ii,jj) * diag_cs%h(ii,jj,k) total_weight = total_weight + weight - ave=ave+field_in(ii,jj,k)*weight + ave = ave+field_in(ii,jj,k)*weight enddo field_out(i,j,k) = ave/(total_weight+eps_face) !Avoid zero mask at all aggregating cells where ave=0.0 - enddo; enddo; enddo - elseif(method == MSK) then !The input field is a mask, subsample + enddo ; enddo ; enddo + elseif (method == MSK) then !The input field is a mask, subsample field_out(:,:,:) = 0.0 - do k= ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d + do k=ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d i0 = isv_o+dl*(i-isv_d) j0 = jsv_o+dl*(j-jsv_d) ave = 0.0 do jj=j0,j0+dl-1 ; do ii=i0,i0+dl-1 - ave=ave+field_in(ii,jj,k) - enddo; enddo - if(ave > 0.0) field_out(i,j,k)=1.0 - enddo; enddo; enddo + ave = ave+field_in(ii,jj,k) + enddo ; enddo + if (ave > 0.0) field_out(i,j,k)=1.0 + enddo ; enddo ; enddo else write (mesg,*) " unknown sampling method: ",method call MOM_error(FATAL, "downsample_field_3d: "//trim(mesg)//" "//trim(diag%debug_str)) @@ -4154,10 +4200,10 @@ subroutine downsample_field_2d(field_in, field_out, dl, method, mask, diag_cs, d weight = mask(ii,jj)*diag_cs%G%areaT(ii,jj) total_weight = total_weight + weight ave = ave+field_in(ii,jj)*weight - enddo; enddo + enddo ; enddo field_out(i,j) = ave/(total_weight + eps_area) !Avoid zero mask at all aggregating cells where ave=0.0 - enddo; enddo - elseif(method == SSP) then ! e.g., T_dfxy_cont_tendency_2d + enddo ; enddo + elseif (method == SSP) then ! e.g., T_dfxy_cont_tendency_2d do j=jsv_d,jev_d ; do i=isv_d,iev_d i0 = isv_o+dl*(i-isv_d) j0 = jsv_o+dl*(j-jsv_d) @@ -4167,11 +4213,11 @@ subroutine downsample_field_2d(field_in, field_out, dl, method, mask, diag_cs, d ! do ii=i0,i0+dl-1 ; do jj=j0,j0+dl-1 weight = mask(ii,jj) total_weight = total_weight + weight - ave=ave+field_in(ii,jj)*weight - enddo; enddo + ave = ave+field_in(ii,jj)*weight + enddo ; enddo field_out(i,j) = ave/(total_weight+epsilon) !Avoid zero mask at all aggregating cells where ave=0.0 - enddo; enddo - elseif(method == PSP) then ! e.g., umo_2d + enddo ; enddo + elseif (method == PSP) then ! e.g., umo_2d do j=jsv_d,jev_d ; do i=isv_d,iev_d i0 = isv_o+dl*(i-isv_d) j0 = jsv_o+dl*(j-jsv_d) @@ -4181,11 +4227,11 @@ subroutine downsample_field_2d(field_in, field_out, dl, method, mask, diag_cs, d do jj=j0,j0+dl-1 weight = mask(ii,jj) total_weight = total_weight +weight - ave=ave+field_in(ii,jj)*weight + ave = ave+field_in(ii,jj)*weight enddo field_out(i,j) = ave/(total_weight+epsilon) !Avoid zero mask at all aggregating cells where ave=0.0 - enddo; enddo - elseif(method == SPP) then ! e.g., vmo_2d + enddo ; enddo + elseif (method == SPP) then ! e.g., vmo_2d do j=jsv_d,jev_d ; do i=isv_d,iev_d i0 = isv_o+dl*(i-isv_d) j0 = jsv_o+dl*(j-jsv_d) @@ -4195,11 +4241,11 @@ subroutine downsample_field_2d(field_in, field_out, dl, method, mask, diag_cs, d do ii=i0,i0+dl-1 weight = mask(ii,jj) total_weight = total_weight +weight - ave=ave+field_in(ii,jj)*weight + ave = ave+field_in(ii,jj)*weight enddo field_out(i,j) = ave/(total_weight+epsilon) !Avoid zero mask at all aggregating cells where ave=0.0 - enddo; enddo - elseif(method == PMP) then + enddo ; enddo + elseif (method == PMP) then do j=jsv_d,jev_d ; do i=isv_d,iev_d i0 = isv_o+dl*(i-isv_d) j0 = jsv_o+dl*(j-jsv_d) @@ -4209,11 +4255,11 @@ subroutine downsample_field_2d(field_in, field_out, dl, method, mask, diag_cs, d do jj=j0,j0+dl-1 weight = mask(ii,jj) * diag_cs%G%dyCu(ii,jj)!*diag_cs%h(ii,jj,1) !Niki? total_weight = total_weight +weight - ave=ave+field_in(ii,jj)*weight + ave = ave+field_in(ii,jj)*weight enddo field_out(i,j) = ave/(total_weight+eps_len) !Avoid zero mask at all aggregating cells where ave=0.0 - enddo; enddo - elseif(method == MPP) then + enddo ; enddo + elseif (method == MPP) then do j=jsv_d,jev_d ; do i=isv_d,iev_d i0 = isv_o+dl*(i-isv_d) j0 = jsv_o+dl*(j-jsv_d) @@ -4223,21 +4269,21 @@ subroutine downsample_field_2d(field_in, field_out, dl, method, mask, diag_cs, d do ii=i0,i0+dl-1 weight = mask(ii,jj)* diag_cs%G%dxCv(ii,jj)!*diag_cs%h(ii,jj,1) !Niki? total_weight = total_weight +weight - ave=ave+field_in(ii,jj)*weight + ave = ave+field_in(ii,jj)*weight enddo field_out(i,j) = ave/(total_weight+eps_len) !Avoid zero mask at all aggregating cells where ave=0.0 - enddo; enddo - elseif(method == MSK) then !The input field is a mask, subsample + enddo ; enddo + elseif (method == MSK) then !The input field is a mask, subsample field_out(:,:) = 0.0 do j=jsv_d,jev_d ; do i=isv_d,iev_d i0 = isv_o+dl*(i-isv_d) j0 = jsv_o+dl*(j-jsv_d) ave = 0.0 do jj=j0,j0+dl-1 ; do ii=i0,i0+dl-1 - ave=ave+field_in(ii,jj) - enddo; enddo - if(ave > 0.0) field_out(i,j)=1.0 - enddo; enddo + ave = ave+field_in(ii,jj) + enddo ; enddo + if (ave > 0.0) field_out(i,j)=1.0 + enddo ; enddo else write (mesg,*) " unknown sampling method: ",method call MOM_error(FATAL, "downsample_field_2d: "//trim(mesg)//" "//trim(diag%debug_str)) @@ -4276,8 +4322,8 @@ subroutine downsample_mask_2d(field_in, field_out, dl, isc_o, jsc_o, isc_d, iec_ do jj=j0,j0+dl-1 ; do ii=i0,i0+dl-1 tot_non_zero = tot_non_zero + field_in(ii,jj) enddo;enddo - if(tot_non_zero > 0.0) field_out(i,j)=1.0 - enddo; enddo + if (tot_non_zero > 0.0) field_out(i,j)=1.0 + enddo ; enddo end subroutine downsample_mask_2d !> Allocate and compute the 3d down sampled mask @@ -4305,15 +4351,15 @@ subroutine downsample_mask_3d(field_in, field_out, dl, isc_o, jsc_o, isc_d, iec_ ks = lbound(field_in,3) ; ke = ubound(field_in,3) allocate(field_out(isd_d:ied_d,jsd_d:jed_d,ks:ke)) field_out(:,:,:) = 0.0 - do k= ks,ke ; do j=jsc_d,jec_d ; do i=isc_d,iec_d + do k=ks,ke ; do j=jsc_d,jec_d ; do i=isc_d,iec_d i0 = isc_o+dl*(i-isc_d) j0 = jsc_o+dl*(j-jsc_d) tot_non_zero = 0.0 do jj=j0,j0+dl-1 ; do ii=i0,i0+dl-1 tot_non_zero = tot_non_zero + field_in(ii,jj,k) enddo;enddo - if(tot_non_zero > 0.0) field_out(i,j,k)=1.0 - enddo; enddo; enddo + if (tot_non_zero > 0.0) field_out(i,j,k)=1.0 + enddo ; enddo ; enddo end subroutine downsample_mask_3d end module MOM_diag_mediator From 772e4563bd78a3c1576cde272342a8e7931ff020 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 30 Nov 2020 06:15:37 -0500 Subject: [PATCH 10/12] Correct units of integrated salt flux diagnostics Corrected the documented units of diagnostics of globally integrated salt flux diagnostics. Also removed two declarations of cmor_units that are redundant with the ordinary units for these diagnostics, and corrected the logical test for an error message related to lateral boundary diffusion. All answers are bitwise identical, but there are changes in the units in available_diags files and in the metadata of some output files. --- src/core/MOM_forcing_type.F90 | 7 +++---- src/parameterizations/lateral/MOM_thickness_diffuse.F90 | 1 - src/parameterizations/vertical/MOM_geothermal.F90 | 2 +- src/tracer/MOM_tracer_hor_diff.F90 | 6 +++--- 4 files changed, 7 insertions(+), 9 deletions(-) diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index ed4b8d1ba2..06f7bf3557 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -1874,17 +1874,16 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, handles%id_total_saltflux = register_scalar_field('ocean_model', & 'total_salt_flux', Time, diag, & - long_name='Area integrated surface salt flux', units='kg', & + long_name='Area integrated surface salt flux', units='kg s-1', & cmor_field_name='total_sfdsi', & - cmor_units='kg s-1', & cmor_standard_name='downward_sea_ice_basal_salt_flux_area_integrated',& cmor_long_name='Downward Sea Ice Basal Salt Flux Area Integrated') handles%id_total_saltFluxIn = register_scalar_field('ocean_model', 'total_salt_Flux_In', & - Time, diag, long_name='Area integrated surface salt flux at surface from coupler', units='kg') + Time, diag, long_name='Area integrated surface salt flux at surface from coupler', units='kg s-1') handles%id_total_saltFluxAdded = register_scalar_field('ocean_model', 'total_salt_Flux_Added', & - Time, diag, long_name='Area integrated surface salt flux due to restoring or flux adjustment', units='kg') + Time, diag, long_name='Area integrated surface salt flux due to restoring or flux adjustment', units='kg s-1') end subroutine register_forcing_type_diags diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index 6ff0184d54..16d9294308 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -2052,7 +2052,6 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS) 'm2 s-1', conversion=US%L_to_m**2*US%s_to_T, & cmor_field_name='diftrblo', & cmor_long_name='Ocean Tracer Diffusivity due to Parameterized Mesoscale Advection', & - cmor_units='m2 s-1', & cmor_standard_name='ocean_tracer_diffusivity_due_to_parameterized_mesoscale_advection') CS%id_KH_u1 = register_diag_field('ocean_model', 'KHTH_u1', diag%axesCu1, Time, & diff --git a/src/parameterizations/vertical/MOM_geothermal.F90 b/src/parameterizations/vertical/MOM_geothermal.F90 index 2ecaa4a78e..dadb5ac526 100644 --- a/src/parameterizations/vertical/MOM_geothermal.F90 +++ b/src/parameterizations/vertical/MOM_geothermal.F90 @@ -581,7 +581,7 @@ subroutine geothermal_init(Time, G, GV, US, param_file, diag, CS, useALEalgorith ! post the static geothermal heating field id = register_static_field('ocean_model', 'geo_heat', diag%axesT1, & 'Geothermal heat flux into ocean', 'W m-2', conversion=US%QRZ_T_to_W_m2, & - cmor_field_name='hfgeou', cmor_units='W m-2', & + cmor_field_name='hfgeou', & cmor_standard_name='upward_geothermal_heat_flux_at_sea_floor', & cmor_long_name='Upward geothermal heat flux at sea floor', & x_cell_method='mean', y_cell_method='mean', area_cell_method='mean') diff --git a/src/tracer/MOM_tracer_hor_diff.F90 b/src/tracer/MOM_tracer_hor_diff.F90 index 43ede7cff5..6506ca4927 100644 --- a/src/tracer/MOM_tracer_hor_diff.F90 +++ b/src/tracer/MOM_tracer_hor_diff.F90 @@ -11,7 +11,7 @@ module MOM_tracer_hor_diff use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type use MOM_domains, only : pass_vector use MOM_debugging, only : hchksum, uvchksum -use MOM_diabatic_driver, only : diabatic_CS +use MOM_diabatic_driver, only : diabatic_CS use MOM_EOS, only : calculate_density, EOS_type, EOS_domain use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe use MOM_error_handler, only : MOM_set_verbosity, callTree_showQuery @@ -60,7 +60,7 @@ module MOM_tracer_hor_diff logical :: use_neutral_diffusion !< If true, use the neutral_diffusion module from within !! tracer_hor_diff. logical :: use_lateral_boundary_diffusion !< If true, use the lateral_boundary_diffusion module from within - !! tracer_hor_diff. + !! tracer_hor_diff. logical :: recalc_neutral_surf !< If true, recalculate the neutral surfaces if CFL has been !! exceeded type(neutral_diffusion_CS), pointer :: neutral_diffusion_CSp => NULL() !< Control structure for neutral diffusion. @@ -1513,7 +1513,7 @@ subroutine tracer_hor_diff_init(Time, G, US, param_file, diag, EOS, diabatic_CSp "USE_NEUTRAL_DIFFUSION and DIFFUSE_ML_TO_INTERIOR are mutually exclusive!") CS%use_lateral_boundary_diffusion = lateral_boundary_diffusion_init(Time, G, param_file, diag, diabatic_CSp, & CS%lateral_boundary_diffusion_CSp) - if (CS%use_neutral_diffusion .and. CS%Diffuse_ML_interior) call MOM_error(FATAL, "MOM_tracer_hor_diff: "// & + if (CS%use_lateral_boundary_diffusion .and. CS%Diffuse_ML_interior) call MOM_error(FATAL, "MOM_tracer_hor_diff: "// & "USE_LATERAL_BOUNDARY_DIFFUSION and DIFFUSE_ML_TO_INTERIOR are mutually exclusive!") call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false.) From d16f04cae2cfdb418978665586fc87d95ddec1da Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Mon, 30 Nov 2020 14:09:34 -0500 Subject: [PATCH 11/12] Adding apt-get update to avoid libpoppler73 issue - Apparently security issues with libpoppler73 have caused the texlive package install to fail. This is fixed by updating the package list before installing any packages. Without texlive-bibtex-extra the `\cite` commands in the doxumentation was causing warnings from doxygen which we consider an error. --- .github/workflows/documentation-and-style.yml | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/.github/workflows/documentation-and-style.yml b/.github/workflows/documentation-and-style.yml index 3988db7675..c83de48159 100644 --- a/.github/workflows/documentation-and-style.yml +++ b/.github/workflows/documentation-and-style.yml @@ -18,11 +18,11 @@ jobs: continue-on-error: true - name: Install packages used when generating documentation - run: > - sudo apt-get install - python3-sphinx python3-lxml perl - texlive-binaries texlive-base bibtool tex-common texlive-bibtex-extra - graphviz + run: | + sudo apt-get update + sudo apt-get install python3-sphinx python3-lxml perl + sudo apt-get install texlive-binaries texlive-base bibtool tex-common texlive-bibtex-extra + sudo apt-get install graphviz - name: Build doxygen HTML run: | From c314fd3d0ff00c039571e19847ca6b01d7d2570a Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 1 Dec 2020 06:28:27 -0500 Subject: [PATCH 12/12] +Always use shorter available_diags format Revised the new available_diags format in response to reviews, including explicitly listing all variable and module names, and always writing the available_diags file in the new shortened format. This includes removing the recently added runtime parameter FULL_AVAILABLE_DIAGS_LIST. All answers are bitwise identical, but there are substantial changes in the available_diags files, and the MOM_parameter_doc files no longer are changed by this PR. --- src/framework/MOM_diag_mediator.F90 | 91 +++++++---------------------- 1 file changed, 20 insertions(+), 71 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index a929613fbc..e04234c859 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -242,9 +242,6 @@ module MOM_diag_mediator !! This file is open if available_diag_doc_unit is > 0. integer :: chksum_iounit = -1 !< The unit number of a diagnostic documentation file. !! This file is open if available_diag_doc_unit is > 0. - logical :: full_diag_list = .true. !< If true, write an entry for every variant of a - !! diagnostic to the available_diags file; otherwise only - !! write an entry for the primary diagnostics. logical :: diag_as_chksum !< If true, log chksums in a text file instead of posting diagnostics logical :: grid_space_axes !< If true, diagnostic horizontal coordinates axes are in grid space. ! The following fields are used for the output of the data. @@ -1989,7 +1986,7 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time integer :: dm_id, i, dl character(len=256) :: msg, cm_string character(len=256) :: new_module_name - character(len=480) :: module_list, var_list, variants + character(len=480) :: module_list, var_list integer :: num_modnm, num_varnm logical :: active @@ -2018,7 +2015,7 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time axes => diag_cs%axesCvi endif - module_list = "{$M" + module_list = "{"//trim(module_name) num_modnm = 1 ! Register the native diagnostic @@ -2033,9 +2030,9 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time y_cell_method=y_cell_method, v_cell_method=v_cell_method, & conversion=conversion, v_extensive=v_extensive) if (associated(axes%xyave_axes)) then - num_varnm = 2 ; var_list = "{$V,$V_xyave" + num_varnm = 2 ; var_list = "{"//trim(field_name)//","//trim(field_name)//"_xyave" else - num_varnm = 1 ; var_list = "{$V" + num_varnm = 1 ; var_list = "{"//trim(field_name) endif if (present(cmor_field_name)) then if (associated(axes%xyave_axes)) then @@ -2090,8 +2087,7 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time if (active) then call diag_remap_set_active(diag_cs%diag_remap_cs(i)) endif - ! module_list = trim(module_list)//","//trim(new_module_name) - module_list = trim(module_list)//",$M_"//trim(diag_cs%diag_remap_cs(i)%diag_module_suffix) + module_list = trim(module_list)//","//trim(new_module_name) num_modnm = num_modnm + 1 endif ! remap_axes%needs_remapping endif ! associated(remap_axes) @@ -2149,8 +2145,7 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time cell_methods=cell_methods, x_cell_method=x_cell_method, & y_cell_method=y_cell_method, v_cell_method=v_cell_method, & conversion=conversion, v_extensive=v_extensive) - ! module_list = trim(module_list)//","//trim(new_module_name) - module_list = trim(module_list)//",$M_2d" + module_list = trim(module_list)//","//trim(new_module_name) num_modnm = num_modnm + 1 endif @@ -2197,8 +2192,7 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time if (active) then call diag_remap_set_active(diag_cs%diag_remap_cs(i)) endif - ! module_list = trim(module_list)//","//trim(new_module_name) - module_list = trim(module_list)//",$M_"//trim(diag_cs%diag_remap_cs(i)%diag_module_suffix)//"_d2" + module_list = trim(module_list)//","//trim(new_module_name) num_modnm = num_modnm + 1 endif ! remap_axes%needs_remapping endif ! associated(remap_axes) @@ -2206,24 +2200,18 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time enddo ! i enddo - if (is_root_pe() .and. (diag_CS%available_diag_doc_unit > 0) .and. .not.diag_CS%full_diag_list) then + if (is_root_pe() .and. (diag_CS%available_diag_doc_unit > 0)) then msg = '' if (present(cmor_field_name)) msg = 'CMOR equivalent is "'//trim(cmor_field_name)//'"' call attach_cell_methods(-1, axes, cm_string, cell_methods, & x_cell_method, y_cell_method, v_cell_method, & v_extensive=v_extensive) module_list = trim(module_list)//"}" - if ((num_modnm > 1) .and. (num_varnm > 1)) then - variants = "mod:"//trim(module_list)//", var:"//trim(var_list) - elseif (num_modnm > 1) then - variants = "mod:"//trim(module_list) - elseif (num_varnm > 1) then - variants = "var:"//trim(var_list) - else - variants = "" - endif - call log_available_diag(dm_id>0, module_name, field_name, cm_string, msg, diag_CS, & - long_name, units, standard_name, variants=variants) + if (num_modnm <= 1) module_list = module_name + if (num_varnm <= 1) var_list = "" + + call log_available_diag(dm_id>0, module_list, field_name, cm_string, msg, diag_CS, & + long_name, units, standard_name, variants=var_list) endif register_diag_field = dm_id @@ -2297,12 +2285,6 @@ logical function register_diag_field_expand_cmor(dm_id, module_name, field_name, call attach_cell_methods(fms_id, axes, cm_string, cell_methods, & x_cell_method, y_cell_method, v_cell_method, & v_extensive=v_extensive) - if (is_root_pe() .and. (diag_CS%available_diag_doc_unit > 0) .and. diag_CS%full_diag_list) then - msg = '' - if (present(cmor_field_name)) msg = 'CMOR equivalent is "'//trim(cmor_field_name)//'"' - call log_available_diag(fms_id>0, module_name, field_name, cm_string, & - msg, diag_CS, long_name, units, standard_name) - endif ! Associated horizontally area-averaged diagnostic fms_xyave_id = DIAG_FIELD_NOT_FOUND if (associated(axes%xyave_axes)) then @@ -2315,12 +2297,6 @@ logical function register_diag_field_expand_cmor(dm_id, module_name, field_name, if (.not. diag_cs%diag_as_chksum) & call attach_cell_methods(fms_xyave_id, axes%xyave_axes, cm_string, & cell_methods, v_cell_method, v_extensive=v_extensive) - if (is_root_pe() .and. (diag_CS%available_diag_doc_unit > 0) .and. diag_CS%full_diag_list) then - msg = '' - if (present(cmor_field_name)) msg = 'CMOR equivalent is "'//trim(cmor_field_name)//'_xyave"' - call log_available_diag(fms_xyave_id>0, module_name, trim(field_name)//'_xyave', cm_string, & - msg, diag_CS, long_name, units, standard_name) - endif endif this_diag => null() if (fms_id /= DIAG_FIELD_NOT_FOUND .or. fms_xyave_id /= DIAG_FIELD_NOT_FOUND) then @@ -2359,12 +2335,6 @@ logical function register_diag_field_expand_cmor(dm_id, module_name, field_name, call attach_cell_methods(fms_id, axes, cm_string, & cell_methods, x_cell_method, y_cell_method, v_cell_method, & v_extensive=v_extensive) - if (is_root_pe() .and. (diag_CS%available_diag_doc_unit > 0) .and. diag_CS%full_diag_list) then - msg = 'native name is "'//trim(field_name)//'"' - call log_available_diag(fms_id>0, module_name, cmor_field_name, cm_string, & - msg, diag_CS, posted_cmor_long_name, posted_cmor_units, & - posted_cmor_standard_name) - endif ! Associated horizontally area-averaged diagnostic fms_xyave_id = DIAG_FIELD_NOT_FOUND if (associated(axes%xyave_axes)) then @@ -2376,12 +2346,6 @@ logical function register_diag_field_expand_cmor(dm_id, module_name, field_name, err_msg=err_msg, interp_method=interp_method, tile_count=tile_count) call attach_cell_methods(fms_xyave_id, axes%xyave_axes, cm_string, & cell_methods, v_cell_method, v_extensive=v_extensive) - if (is_root_pe() .and. (diag_CS%available_diag_doc_unit > 0) .and. diag_CS%full_diag_list) then - msg = 'native name is "'//trim(field_name)//'_xyave"' - call log_available_diag(fms_xyave_id>0, module_name, trim(cmor_field_name)//'_xyave', & - cm_string, msg, diag_CS, posted_cmor_long_name, posted_cmor_units, & - posted_cmor_standard_name) - endif endif this_diag => null() if (fms_id /= DIAG_FIELD_NOT_FOUND .or. fms_xyave_id /= DIAG_FIELD_NOT_FOUND) then @@ -2797,19 +2761,14 @@ function register_scalar_field(module_name, field_name, init_time, diag_cs, & ! Document diagnostics in list of available diagnostics if (is_root_pe() .and. diag_CS%available_diag_doc_unit > 0) then - if (present(cmor_field_name) .and. .not.diag_CS%full_diag_list) then + if (present(cmor_field_name)) then call log_available_diag(associated(diag), module_name, field_name, '', '', diag_CS, & long_name, units, standard_name, & - variants="var:{$V,"//trim(cmor_field_name)//"}") + variants="{"//trim(field_name)//","//trim(cmor_field_name)//"}") else call log_available_diag(associated(diag), module_name, field_name, '', '', diag_CS, & long_name, units, standard_name) endif - if (present(cmor_field_name) .and. diag_CS%full_diag_list) then - call log_available_diag(associated(cmor_diag), module_name, cmor_field_name, & - '', '', diag_CS, posted_cmor_long_name, posted_cmor_units, & - posted_cmor_standard_name) - endif endif register_scalar_field = dm_id @@ -2950,19 +2909,14 @@ function register_static_field(module_name, field_name, axes, & ! Document diagnostics in list of available diagnostics if (is_root_pe() .and. diag_CS%available_diag_doc_unit > 0) then - if (present(cmor_field_name) .and. .not.diag_CS%full_diag_list) then + if (present(cmor_field_name)) then call log_available_diag(associated(diag), module_name, field_name, '', '', diag_CS, & long_name, units, standard_name, & - variants="var:{$V,"//trim(cmor_field_name)//"}") + variants="{"//trim(field_name)//","//trim(cmor_field_name)//"}") else call log_available_diag(associated(diag), module_name, field_name, '', '', diag_CS, & long_name, units, standard_name) endif - if (present(cmor_field_name) .and. diag_CS%full_diag_list) then - call log_available_diag(associated(cmor_diag), module_name, cmor_field_name, & - '', '', diag_CS, posted_cmor_long_name, posted_cmor_units, & - posted_cmor_standard_name) - endif endif register_static_field = dm_id @@ -3226,12 +3180,6 @@ subroutine diag_mediator_init(G, GV, US, nz, param_file, diag_cs, doc_file_dir) if ((.not.opened) .or. (ios /= 0)) then call MOM_error(FATAL, "Failed to open available diags file "//trim(doc_path)//".") endif - - call get_param(param_file, mdl, "FULL_AVAILABLE_DIAGS_LIST", diag_CS%full_diag_list, & - "If true, write entries for every variant of a diagnostic to the available //"& - "diags file; otherwise only write an entry for the primary diagnostics."//& - "ocean diagnostics that can be included in a diag_table.", & - default=.true.) endif endif @@ -3598,15 +3546,16 @@ subroutine log_available_diag(used, module_name, field_name, cell_methods_string character(len=240) :: mesg if (used) then - mesg = '"'//trim(module_name)//'", "'//trim(field_name)//'" [Used]' + mesg = '"'//trim(field_name)//'" [Used]' else - mesg = '"'//trim(module_name)//'", "'//trim(field_name)//'" [Unused]' + mesg = '"'//trim(field_name)//'" [Unused]' endif if (len(trim((comment)))>0) then write(diag_CS%available_diag_doc_unit, '(a,x,"(",a,")")') trim(mesg),trim(comment) else write(diag_CS%available_diag_doc_unit, '(a)') trim(mesg) endif + call describe_option("modules", module_name, diag_CS) if (present(long_name)) call describe_option("long_name", long_name, diag_CS) if (present(units)) call describe_option("units", units, diag_CS) if (present(standard_name)) &