From dc7476755e65e8d594f2cc8913c6b8b23b69a80b Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Mon, 11 Jan 2021 21:58:42 -0500 Subject: [PATCH 01/11] Adds targets check_mom6_api_nuopc, check_mom6_api_coupled - Allows us to check that we have not broken the ability to compile with the NUOPC, MCT or coupled_driver modules. - This is not a complete compile but only upto the last module that does not need anything other than MOM6 objects. --- .testing/Makefile | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/.testing/Makefile b/.testing/Makefile index 33d3ea4c4e..a780f19323 100644 --- a/.testing/Makefile +++ b/.testing/Makefile @@ -328,6 +328,40 @@ $(DEPS)/Makefile: ../ac/deps/Makefile mkdir -p $(@D) cp $< $@ +#--- +# The following block does a non-library build of a coupled driver interface to MOM, along with everything below it. +# This simply checks that we have not broken the ability to compile. This is not a means to build a complete coupled executable. +# Todo: +# - avoid re-building FMS and MOM6 src by re-using existing object/mod files +# - use autoconf rather than mkmf templates +MK_TEMPLATE ?= ../../$(DEPS)/mkmf/templates/ncrc-gnu.mk +# NUOPC driver +build/nuopc/Makefile: ../config_src/nuopc_driver $(MOM_SOURCE) + mkdir -p $(@D) + cd $(@D) \ + && $(MOM_ENV) ../../$(DEPS)/mkmf/bin/list_paths -l ../../$(DEPS)/fms/src ../../../config_src/nuopc_driver ../../../config_src/dynamic_symmetric ../../../src/ ../../../config_src/external \ + && $(MOM_ENV) ../../$(DEPS)/mkmf/bin/mkmf -t $(MK_TEMPLATE) -p MOM6 path_names +build/nuopc/mom_ocean_model_nuopc.o: build/nuopc/Makefile + cd $(@D) && make $(@F) +check_mom6_api_nuopc: build/nuopc/mom_ocean_model_nuopc.o +# GFDL coupled driver +build/coupled/Makefile: ../config_src/coupled_driver $(MOM_SOURCE) + mkdir -p $(@D) + cd $(@D) \ + && $(MOM_ENV) ../../$(DEPS)/mkmf/bin/list_paths -l ../../$(DEPS)/fms/src ../../../config_src/coupled_driver ../../../config_src/dynamic_symmetric ../../../src/ ../../../config_src/external \ + && $(MOM_ENV) ../../$(DEPS)/mkmf/bin/mkmf -t $(MK_TEMPLATE) -p MOM6 path_names +build/coupled/ocean_model_MOM.o: build/coupled/Makefile + cd $(@D) && make $(@F) +check_mom6_api_coupled: build/coupled/ocean_model_MOM.o +# MCT driver +build/mct/Makefile: ../config_src/mct_driver $(MOM_SOURCE) + mkdir -p $(@D) + cd $(@D) \ + && $(MOM_ENV) ../../$(DEPS)/mkmf/bin/list_paths -l ../../$(DEPS)/fms/src ../../../config_src/mct_driver ../../../config_src/dynamic_symmetric ../../../src/ ../../../config_src/external \ + && $(MOM_ENV) ../../$(DEPS)/mkmf/bin/mkmf -t $(MK_TEMPLATE) -p MOM6 path_names +build/mct/mom_ocean_model_mct.o: build/mct/Makefile + cd $(@D) && make $(@F) +check_mom6_api_mct: build/mct/mom_ocean_model_mct.o #--- # Python preprocessing From 5c93def092eae6f06df4cd94bc44de23550cf36d Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Tue, 12 Jan 2021 14:14:21 -0500 Subject: [PATCH 02/11] MOM_hor_visc: horizontal_viscosity loop reorder This patch reorders many of the loops in horizontal_viscosity in order to improve vectorization of the Laplacian and biharmonic viscosities. Specifically, a single loop containing many different computations were broken up into many loops of individual operations. This patch required introduction of several new 2D arrays. On Gaea's Broadwell CPUs (E5-2697 v4), this is a ~80% speedup on a 32x32x75 `benchmark` configuration. Smaller speedups were observed on AMD processors. On the Gaea nodes, performance appears to be limited by the very large number of variables in the function stack, and the high degree of stack spill. Further loop reordering may cause slowdowns unless the stack usage is reduced. No answers should be changed by this patch, but deserves extra scrutiny given the fundamental role of this function in nearly all simulations. --- .../lateral/MOM_hor_visc.F90 | 689 ++++++++++++------ 1 file changed, 463 insertions(+), 226 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index ed3ef7173e..06a6269dc5 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -196,7 +196,6 @@ module MOM_hor_visc integer :: id_normstress = -1, id_shearstress = -1 !>@} - end type hor_visc_CS contains @@ -265,8 +264,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, str_xx_GME,& ! smoothed diagonal term in the stress tensor from GME [H L2 T-2 ~> m3 s-2 or kg s-2] bhstr_xx, & ! A copy of str_xx that only contains the biharmonic contribution [H L2 T-2 ~> m3 s-2 or kg s-2] FrictWorkIntz, & ! depth integrated energy dissipated by lateral friction [R L2 T-3 ~> W m-2] - ! Leith_Kh_h, & ! Leith Laplacian viscosity at h-points [L2 T-1 ~> m2 s-1] - ! Leith_Ah_h, & ! Leith bi-harmonic viscosity at h-points [L4 T-1 ~> m4 s-1] grad_vort_mag_h, & ! Magnitude of vorticity gradient at h-points [L-1 T-1 ~> m-1 s-1] grad_vort_mag_h_2d, & ! Magnitude of 2d vorticity gradient at h-points [L-1 T-1 ~> m-1 s-1] Del2vort_h, & ! Laplacian of vorticity at h-points [L-2 T-1 ~> m-2 s-1] @@ -289,7 +286,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, str_xy, & ! str_xy is the cross term in the stress tensor [H L2 T-2 ~> m3 s-2 or kg s-2] str_xy_GME, & ! smoothed cross term in the stress tensor from GME [H L2 T-2 ~> m3 s-2 or kg s-2] bhstr_xy, & ! A copy of str_xy that only contains the biharmonic contribution [H L2 T-2 ~> m3 s-2 or kg s-2] - vort_xy, & ! Vertical vorticity (dv/dx - du/dy) including metric terms [T-1 ~> s-1] + vort_xy, & ! Vertical vorticity (dv/dx - du/dy) including metric terms [T-1 ~> s-1] Leith_Kh_q, & ! Leith Laplacian viscosity at q-points [L2 T-1 ~> m2 s-1] Leith_Ah_q, & ! Leith bi-harmonic viscosity at q-points [L4 T-1 ~> m4 s-1] grad_vort_mag_q, & ! Magnitude of vorticity gradient at q-points [L-1 T-1 ~> m-1 s-1] @@ -297,8 +294,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, Del2vort_q, & ! Laplacian of vorticity at q-points [L-2 T-1 ~> m-2 s-1] grad_div_mag_q, & ! Magnitude of divergence gradient at q-points [L-1 T-1 ~> m-1 s-1] grad_vel_mag_q, & ! Magnitude of the velocity gradient tensor squared at q-points [T-2 ~> s-2] - hq, & ! harmonic mean of the harmonic means of the u- & v point thicknesses [H ~> m or kg m-2] - ! This form guarantees that hq/hu < 4. + hq, & ! harmonic mean of the harmonic means of the u- & v point thicknesses [H ~> m or kg m-2] + ! This form guarantees that hq/hu < 4. grad_vel_mag_bt_q, & ! Magnitude of the barotropic velocity gradient tensor squared at q-points [T-2 ~> s-2] boundary_mask_q ! A mask that zeroes out cells with at least one land edge [nondim] @@ -323,30 +320,25 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, div_xx_h, & ! horizontal divergence [T-1 ~> s-1] sh_xx_h, & ! horizontal tension (du/dx - dv/dy) including metric terms [T-1 ~> s-1] NoSt ! A diagnostic array of normal stress [T-1 ~> s-1]. - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & - grid_Re_Kh, & !< Grid Reynolds number for Laplacian horizontal viscosity at h points [nondim] - grid_Re_Ah, & !< Grid Reynolds number for Biharmonic horizontal viscosity at h points [nondim] - GME_coeff_h !< GME coeff. at h-points [L2 T-1 ~> m2 s-1] - real :: Ah ! biharmonic viscosity [L4 T-1 ~> m4 s-1] - real :: Kh ! Laplacian viscosity [L2 T-1 ~> m2 s-1] + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & + grid_Re_Kh, & ! Grid Reynolds number for Laplacian horizontal viscosity at h points [nondim] + grid_Re_Ah, & ! Grid Reynolds number for Biharmonic horizontal viscosity at h points [nondim] + GME_coeff_h ! GME coeff. at h-points [L2 T-1 ~> m2 s-1] real :: AhSm ! Smagorinsky biharmonic viscosity [L4 T-1 ~> m4 s-1] real :: AhLth ! 2D Leith biharmonic viscosity [L4 T-1 ~> m4 s-1] real :: mod_Leith ! nondimensional coefficient for divergence part of modified Leith ! viscosity. Here set equal to nondimensional Laplacian Leith constant. ! This is set equal to zero if modified Leith is not used. - real :: Shear_mag ! magnitude of the shear [T-1 ~> s-1] - real :: vert_vort_mag ! magnitude of the vertical vorticity gradient [L-1 T-1 ~> m-1 s-1] + real :: Shear_mag_bc ! Shear_mag value in backscatter [T-1 ~> s-1] + real :: sh_xx_sq ! Square of tension (sh_xx) [T-2 ~> s-2] + real :: sh_xy_sq ! Square of shearing strain (sh_xy) [T-2 ~> s-2] real :: h2uq, h2vq ! temporary variables [H2 ~> m2 or kg2 m-4]. real :: hu, hv ! Thicknesses interpolated by arithmetic means to corner ! points; these are first interpolated to u or v velocity ! points where masks are applied [H ~> m or kg m-2]. real :: h_neglect ! thickness so small it can be lost in roundoff and so neglected [H ~> m or kg m-2] real :: h_neglect3 ! h_neglect^3 [H3 ~> m3 or kg3 m-6] - real :: hrat_min ! minimum thicknesses at the 4 neighboring - ! velocity points divided by the thickness at the stress - ! point (h or q point) [nondim] - real :: visc_bound_rem ! fraction of overall viscous bounds that - ! remain to be applied [nondim] + real :: h_min ! Minimum h at the 4 neighboring velocity points [H ~> m] real :: Kh_scale ! A factor between 0 and 1 by which the horizontal ! Laplacian viscosity is rescaled [nondim] real :: RoScl ! The scaling function for MEKE source term [nondim] @@ -365,6 +357,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! calculation gives the same value as if f were 0 [nondim]. real :: H0_GME ! Depth used to scale down GME coefficient in shallow areas [Z ~> m] real :: KE ! Local kinetic energy [L2 T-2 ~> m2 s-2] + real :: d_del2u ! dy-weighted Laplacian(u) diff in x [L-2 T-1 ~> m-2 s-1] + real :: d_del2v ! dx-weighted Laplacian(v) diff in y [L-2 T-1 ~> m-2 s-1] + real :: d_str ! Stress tensor update [H L2 T-2 ~> m3 s-2 or kg s-2] + real :: grad_vort ! Vorticity gradient magnitude [L-1 T-1 ~> m-1 s-1] + real :: grad_vort_qg ! QG-based vorticity gradient magnitude [L-1 T-1 ~> m-1 s-1] + real :: grid_Kh ! Laplacian viscosity bound by grid [L2 T-1 ~> m2 s-1] + real :: grid_Ah ! Biharmonic viscosity bound by grid [L4 T-1 ~> m4 s-1] logical :: rescale_Kh, legacy_bound logical :: find_FrictWork @@ -374,6 +373,18 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz integer :: i, j, k, n real :: inv_PI3, inv_PI2, inv_PI6 + + ! Fields evaluated on active layers, used for constructing 3D stress fields + ! NOTE: The position of these declarations can impact performance, due to the + ! very large number of stack arrays in this function. Move with caution! + real, dimension(SZIB_(G),SZJB_(G)) :: & + Ah, & ! biharmonic viscosity (h or q) [L4 T-1 ~> m4 s-1] + Kh, & ! Laplacian viscosity [L2 T-1 ~> m2 s-1] + Shear_mag, & ! magnitude of the shear [T-1 ~> s-1] + vert_vort_mag, & ! magnitude of the vertical vorticity gradient [L-1 T-1 ~> m-1 s-1] + hrat_min, & ! h_min divided by the thickness at the stress point (h or q) [nondim] + visc_bound_rem ! fraction of overall viscous bounds that remain to be applied [nondim] + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB @@ -383,9 +394,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, inv_PI2 = 1.0/((4.0*atan(1.0))**2) inv_PI6 = inv_PI3 * inv_PI3 - Ah_h(:,:,:) = 0.0 - Kh_h(:,:,:) = 0.0 - if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%OBC_pe) then apply_OBC = OBC%Flather_u_BCs_exist_globally .or. OBC%Flather_v_BCs_exist_globally apply_OBC = .true. @@ -505,10 +513,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, !$OMP str_xx, str_xy, bhstr_xx, bhstr_xy, str_xx_GME, str_xy_GME, & !$OMP vort_xy, vort_xy_dx, vort_xy_dy, div_xx, div_xx_dx, div_xx_dy, & !$OMP grad_div_mag_h, grad_div_mag_q, grad_vort_mag_h, grad_vort_mag_q, & - !$OMP grad_vort_mag_h_2d, grad_vort_mag_q_2d, & + !$OMP grad_vort, grad_vort_qg, grad_vort_mag_h_2d, grad_vort_mag_q_2d, & !$OMP grad_vel_mag_h, grad_vel_mag_q, & !$OMP grad_vel_mag_bt_h, grad_vel_mag_bt_q, grad_d2vel_mag_h, & - !$OMP meke_res_fn, Shear_mag, vert_vort_mag, hrat_min, visc_bound_rem, & + !$OMP meke_res_fn, Shear_mag, Shear_mag_bc, vert_vort_mag, h_min, hrat_min, visc_bound_rem, & + !$OMP sh_xx_sq, sh_xy_sq, grid_Ah, grid_Kh, d_Del2u, d_Del2v, d_str, & !$OMP Kh, Ah, AhSm, AhLth, local_strain, Sh_F_pow, & !$OMP dDel2vdx, dDel2udy, DY_dxCv, DX_dyCu, Del2vort_q, Del2vort_h, KE, & !$OMP h2uq, h2vq, hu, hv, hq, FatH, RoScl, GME_coeff & @@ -519,22 +528,25 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! shearing strain advocated by Smagorinsky (1993) and discussed in ! Griffies and Hallberg (2000). - ! Calculate horizontal tension - do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2 + do j=Jsq-2,Jeq+2 ; do i=Isq-2,Ieq+2 + ! Calculate horizontal tension dudx(i,j) = CS%DY_dxT(i,j)*(G%IdyCu(I,j) * u(I,j,k) - & G%IdyCu(I-1,j) * u(I-1,j,k)) dvdy(i,j) = CS%DX_dyT(i,j)*(G%IdxCv(i,J) * v(i,J,k) - & G%IdxCv(i,J-1) * v(i,J-1,k)) sh_xx(i,j) = dudx(i,j) - dvdy(i,j) - if (CS%id_normstress > 0) NoSt(i,j,k) = sh_xx(i,j) - enddo ; enddo - ! Components for the shearing strain - do J=Jsq-2,Jeq+2 ; do I=Isq-2,Ieq+2 + ! Components for the shearing strain dvdx(I,J) = CS%DY_dxBu(I,J)*(v(i+1,J,k)*G%IdyCv(i+1,J) - v(i,J,k)*G%IdyCv(i,J)) dudy(I,J) = CS%DX_dyBu(I,J)*(u(I,j+1,k)*G%IdxCu(I,j+1) - u(I,j,k)*G%IdxCu(I,j)) enddo ; enddo + if (CS%id_normstress > 0) then + do j=Jsq-2,Jeq+2 ; do i=Isq-2,Ieq+2 + NoSt(i,j,k) = sh_xx(i,j) + enddo ; enddo + endif + ! Interpolate the thicknesses to velocity points. ! The extra wide halos are to accommodate the cross-corner-point projections ! in OBCs, which are not ordinarily be necessary, and might not be necessary @@ -608,7 +620,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, endif endif - if (OBC%segment(n)%direction == OBC_DIRECTION_N) then ! There are extra wide halos here to accommodate the cross-corner-point ! OBC projections, but they might not be necessary if the accelerations @@ -765,7 +776,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, grad_div_mag_h(i,j) =sqrt((0.5*(div_xx_dx(I,j) + div_xx_dx(I-1,j)))**2 + & (0.5 * (div_xx_dy(i,J) + div_xx_dy(i,J-1)))**2) enddo ; enddo - !do J=js-1,Jeq ; do I=is-1,Ieq do j=Jsq-1,Jeq+1 ; do i=Isq-1,Ieq+1 grad_div_mag_q(I,J) =sqrt((0.5*(div_xx_dx(I,j) + div_xx_dx(I,j+1)))**2 + & (0.5 * (div_xx_dy(i,J) + div_xx_dy(i+1,J)))**2) @@ -828,133 +838,247 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, meke_res_fn = 1. - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - if ((CS%Smagorinsky_Kh) .or. (CS%Smagorinsky_Ah)) then - Shear_mag = sqrt(sh_xx(i,j)*sh_xx(i,j) + & - 0.25*((sh_xy(I-1,J-1)*sh_xy(I-1,J-1) + sh_xy(I,J)*sh_xy(I,J)) + & - (sh_xy(I-1,J)*sh_xy(I-1,J) + sh_xy(I,J-1)*sh_xy(I,J-1)))) + if ((CS%Smagorinsky_Kh) .or. (CS%Smagorinsky_Ah)) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + sh_xx_sq = sh_xx(i,j) * sh_xx(i,j) + sh_xy_sq = 0.25 * ( & + (sh_xy(I-1,J-1) * sh_xy(I-1,J-1) + sh_xy(I,J) * sh_xy(I,J)) & + + (sh_xy(I-1,J) * sh_xy(I-1,J) + sh_xy(I,J-1) * sh_xy(I,J-1)) & + ) + Shear_mag(i,j) = sqrt(sh_xx_sq + sh_xy_sq) + enddo ; enddo + endif + + if (CS%better_bound_Ah .or. CS%better_bound_Kh) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + h_min = min(h_u(I,j), h_u(I-1,j), h_v(i,J), h_v(i,J-1)) + hrat_min(i,j) = min(1.0, h_min / (h(i,j,k) + h_neglect)) + enddo ; enddo + + if (CS%better_bound_Kh) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + visc_bound_rem(i,j) = 1.0 + enddo ; enddo endif + endif + + if (CS%Laplacian) then if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then if (CS%use_QG_Leith_visc) then - vert_vort_mag = MIN(grad_vort_mag_h(i,j) + grad_div_mag_h(i,j),3.*grad_vort_mag_h_2d(i,j)) + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + grad_vort = grad_vort_mag_h(i,j) + grad_div_mag_h(i,j) + grad_vort_qg = 3. * grad_vort_mag_h_2d(i,j) + vert_vort_mag(i,j) = min(grad_vort, grad_vort_qg) + enddo ; enddo else - vert_vort_mag = (grad_vort_mag_h(i,j) + grad_div_mag_h(i,j)) + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + vert_vort_mag(i,j) = grad_vort_mag_h(i,j) + grad_div_mag_h(i,j) + enddo ; enddo endif endif - if (CS%better_bound_Ah .or. CS%better_bound_Kh) then - hrat_min = min(1.0, min(h_u(I,j), h_u(I-1,j), h_v(i,J), h_v(i,J-1)) / & - (h(i,j,k) + h_neglect) ) - visc_bound_rem = 1.0 - endif - if (CS%Laplacian) then - ! Determine the Laplacian viscosity at h points, using the - ! largest value from several parameterizations. - Kh = CS%Kh_bg_xx(i,j) ! Static (pre-computed) background viscosity + ! Determine the Laplacian viscosity at h points, using the + ! largest value from several parameterizations. + + ! Static (pre-computed) background viscosity + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Kh(i,j) = CS%Kh_bg_xx(i,j) + enddo ; enddo + + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 if (CS%add_LES_viscosity) then - if (CS%Smagorinsky_Kh) Kh = Kh + CS%Laplac2_const_xx(i,j) * Shear_mag - if (CS%Leith_Kh) Kh = Kh + CS%Laplac3_const_xx(i,j) * vert_vort_mag*inv_PI3 + if (CS%Smagorinsky_Kh) & + Kh(i,j) = Kh(i,j) + CS%Laplac2_const_xx(i,j) * Shear_mag(i,j) + if (CS%Leith_Kh) & + Kh(i,j) = Kh(i,j) + CS%Laplac3_const_xx(i,j) * vert_vort_mag(i,j) * inv_PI3 else - if (CS%Smagorinsky_Kh) Kh = max( Kh, CS%Laplac2_const_xx(i,j) * Shear_mag ) - if (CS%Leith_Kh) Kh = max( Kh, CS%Laplac3_const_xx(i,j) * vert_vort_mag*inv_PI3) + if (CS%Smagorinsky_Kh) & + Kh(i,j) = max(Kh(i,j), CS%Laplac2_const_xx(i,j) * Shear_mag(i,j)) + if (CS%Leith_Kh) & + Kh(i,j) = max(Kh(i,j), CS%Laplac3_const_xx(i,j) * vert_vort_mag(i,j) * inv_PI3) endif - ! All viscosity contributions above are subject to resolution scaling - if (rescale_Kh) Kh = VarMix%Res_fn_h(i,j) * Kh - if (CS%res_scale_MEKE) meke_res_fn = VarMix%Res_fn_h(i,j) + enddo ; enddo + + ! All viscosity contributions above are subject to resolution scaling + + if (rescale_Kh) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Kh(i,j) = VarMix%Res_fn_h(i,j) * Kh(i,j) + enddo ; enddo + endif + + if (legacy_bound) then ! Older method of bounding for stability - if (legacy_bound) Kh = min(Kh, CS%Kh_Max_xx(i,j)) - Kh = max( Kh, CS%Kh_bg_min ) ! Place a floor on the viscosity, if desired. + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Kh(i,j) = min(Kh(i,j), CS%Kh_Max_xx(i,j)) + enddo ; enddo + endif + + ! Place a floor on the viscosity, if desired. + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Kh(i,j) = max(Kh(i,j), CS%Kh_bg_min) + enddo ; enddo + + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + if (CS%res_scale_MEKE) meke_res_fn = VarMix%Res_fn_h(i,j) + if (use_MEKE_Ku) & - Kh = Kh + MEKE%Ku(i,j) * meke_res_fn ! *Add* the MEKE contribution (might be negative) - if (CS%anisotropic) Kh = Kh + CS%Kh_aniso * ( 1. - CS%n1n2_h(i,j)**2 ) ! *Add* the tension component - ! of anisotropic viscosity + ! *Add* the MEKE contribution (might be negative) + Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j) * meke_res_fn + enddo ; enddo - ! Newer method of bounding for stability + if (CS%anisotropic) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + ! *Add* the tension component of anisotropic viscosity + Kh(i,j) = Kh(i,j) + CS%Kh_aniso * (1. - CS%n1n2_h(i,j)**2) + enddo ; enddo + endif + + ! Newer method of bounding for stability + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 if (CS%better_bound_Kh) then - if (Kh >= hrat_min*CS%Kh_Max_xx(i,j)) then - visc_bound_rem = 0.0 - Kh = hrat_min*CS%Kh_Max_xx(i,j) + if (Kh(i,j) >= hrat_min(i,j) * CS%Kh_Max_xx(i,j)) then + visc_bound_rem(i,j) = 0.0 + Kh(i,j) = hrat_min(i,j) * CS%Kh_Max_xx(i,j) else - visc_bound_rem = 1.0 - Kh / (hrat_min*CS%Kh_Max_xx(i,j)) + visc_bound_rem(i,j) = 1.0 - Kh(i,j) / (hrat_min(i,j) * CS%Kh_Max_xx(i,j)) endif endif + enddo ; enddo - if ((CS%id_Kh_h>0) .or. find_FrictWork .or. CS%debug) Kh_h(i,j,k) = Kh + if (CS%id_Kh_h>0 .or. CS%debug) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Kh_h(i,j,k) = Kh(i,j) + enddo ; enddo + endif - if (CS%id_grid_Re_Kh>0) then + if (CS%id_grid_Re_Kh>0) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 KE = 0.125*((u(I,j,k)+u(I-1,j,k))**2 + (v(i,J,k)+v(i,J-1,k))**2) - grid_Re_Kh(i,j,k) = (sqrt(KE) * sqrt(CS%grid_sp_h2(i,j))) & - / max(Kh, CS%min_grid_Kh) - endif + grid_Kh = max(Kh(i,j), CS%min_grid_Kh) + grid_Re_Kh(i,j,k) = (sqrt(KE) * sqrt(CS%grid_sp_h2(i,j))) / grid_Kh + enddo ; enddo + endif - if (CS%id_div_xx_h>0) div_xx_h(i,j,k) = div_xx(i,j) - if (CS%id_sh_xx_h>0) sh_xx_h(i,j,k) = sh_xx(i,j) + if (CS%id_div_xx_h>0) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + div_xx_h(i,j,k) = div_xx(i,j) + enddo ; enddo + endif + + if (CS%id_sh_xx_h>0) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + sh_xx_h(i,j,k) = sh_xx(i,j) + enddo ; enddo + endif - str_xx(i,j) = -Kh * sh_xx(i,j) - else ! not Laplacian + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + str_xx(i,j) = -Kh(i,j) * sh_xx(i,j) + enddo ; enddo + else + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 str_xx(i,j) = 0.0 - endif ! Laplacian + enddo ; enddo + endif - if (CS%anisotropic) then + if (CS%anisotropic) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 ! Shearing-strain averaged to h-points local_strain = 0.25 * ( (sh_xy(I,J) + sh_xy(I-1,J-1)) + (sh_xy(I-1,J) + sh_xy(I,J-1)) ) ! *Add* the shear-strain contribution to the xx-component of stress str_xx(i,j) = str_xx(i,j) - CS%Kh_aniso * CS%n1n2_h(i,j) * CS%n1n1_m_n2n2_h(i,j) * local_strain - endif + enddo ; enddo + endif - if (CS%biharmonic) then - ! Determine the biharmonic viscosity at h points, using the - ! largest value from several parameterizations. - AhSm = 0.0; AhLth = 0.0 - if ((CS%Smagorinsky_Ah) .or. (CS%Leith_Ah)) then - if (CS%Smagorinsky_Ah) then - if (CS%bound_Coriolis) then - AhSm = Shear_mag * (CS%Biharm_const_xx(i,j) + & - CS%Biharm_const2_xx(i,j)*Shear_mag) - else - AhSm = CS%Biharm_const_xx(i,j) * Shear_mag - endif - endif - if (CS%Leith_Ah) AhLth = CS%Biharm6_const_xx(i,j) * abs(Del2vort_h(i,j)) * inv_PI6 - Ah = MAX(MAX(CS%Ah_bg_xx(i,j), AhSm), AhLth) - if (CS%bound_Ah .and. .not.CS%better_bound_Ah) & - Ah = MIN(Ah, CS%Ah_Max_xx(i,j)) - else - Ah = CS%Ah_bg_xx(i,j) - endif ! Smagorinsky_Ah or Leith_Ah + if (CS%biharmonic) then + ! Determine the biharmonic viscosity at h points, using the + ! largest value from several parameterizations. + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Ah(i,j) = CS%Ah_bg_xx(i,j) + enddo ; enddo - if (use_MEKE_Au) Ah = Ah + MEKE%Au(i,j) ! *Add* the MEKE contribution + if ((CS%Smagorinsky_Ah) .or. (CS%Leith_Ah)) then + if (CS%Smagorinsky_Ah) then + if (CS%bound_Coriolis) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + AhSm = Shear_mag(i,j) * (CS%Biharm_const_xx(i,j) & + + CS%Biharm_const2_xx(i,j) * Shear_mag(i,j) & + ) + Ah(i,j) = max(Ah(i,j), AhSm) + enddo ; enddo + else + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + AhSm = CS%Biharm_const_xx(i,j) * Shear_mag(i,j) + Ah(i,j) = max(Ah(i,j), AhSm) + enddo ; enddo + endif + endif - if (CS%Re_Ah > 0.0) then - KE = 0.125*((u(I,j,k)+u(I-1,j,k))**2 + (v(i,J,k)+v(i,J-1,k))**2) - Ah = sqrt(KE) * CS%Re_Ah_const_xx(i,j) + if (CS%Leith_Ah) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + AhLth = CS%Biharm6_const_xx(i,j) * abs(Del2vort_h(i,j)) * inv_PI6 + Ah(i,j) = max(Ah(i,j), AhLth) + enddo ; enddo endif - if (CS%better_bound_Ah) then - Ah = MIN(Ah, visc_bound_rem*hrat_min*CS%Ah_Max_xx(i,j)) + if (CS%bound_Ah .and. .not. CS%better_bound_Ah) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Ah(i,j) = min(Ah(i,j), CS%Ah_Max_xx(i,j)) + enddo ; enddo endif + endif ! Smagorinsky_Ah or Leith_Ah - if ((CS%id_Ah_h>0) .or. find_FrictWork .or. CS%debug) Ah_h(i,j,k) = Ah + if (use_MEKE_Au) then + ! *Add* the MEKE contribution + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Ah(i,j) = Ah(i,j) + MEKE%Au(i,j) + enddo ; enddo + endif - if (CS%id_grid_Re_Ah>0) then + if (CS%Re_Ah > 0.0) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 KE = 0.125*((u(I,j,k)+u(I-1,j,k))**2 + (v(i,J,k)+v(i,J-1,k))**2) - grid_Re_Ah(i,j,k) = (sqrt(KE) * CS%grid_sp_h3(i,j)) & - / max(Ah, CS%min_grid_Ah) + Ah(i,j) = sqrt(KE) * CS%Re_Ah_const_xx(i,j) + enddo ; enddo + endif + + if (CS%better_bound_Ah) then + if (CS%better_bound_Kh) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Ah(i,j) = min(Ah(i,j), visc_bound_rem(i,j) * hrat_min(i,j) * CS%Ah_Max_xx(i,j)) + enddo ; enddo + else + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Ah(i,j) = min(Ah(i,j), hrat_min(i,j) * CS%Ah_Max_xx(i,j)) + enddo ; enddo endif + endif - str_xx(i,j) = str_xx(i,j) + Ah * & - (CS%DY_dxT(i,j) * (G%IdyCu(I,j)*Del2u(I,j) - G%IdyCu(I-1,j)*Del2u(I-1,j)) - & - CS%DX_dyT(i,j) * (G%IdxCv(i,J)*Del2v(i,J) - G%IdxCv(i,J-1)*Del2v(i,J-1))) + if ((CS%id_Ah_h>0) .or. CS%debug) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Ah_h(i,j,k) = Ah(i,j) + enddo ; enddo + endif - ! Keep a copy of the biharmonic contribution for backscatter parameterization - bhstr_xx(i,j) = Ah * & - (CS%DY_dxT(i,j) * (G%IdyCu(I,j)*Del2u(I,j) - G%IdyCu(I-1,j)*Del2u(I-1,j)) - & - CS%DX_dyT(i,j) * (G%IdxCv(i,J)*Del2v(i,J) - G%IdxCv(i,J-1)*Del2v(i,J-1))) - bhstr_xx(i,j) = bhstr_xx(i,j) * (h(i,j,k) * CS%reduction_xx(i,j)) + if (CS%id_grid_Re_Ah>0) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + KE = 0.125 * ((u(I,j,k) + u(I-1,j,k))**2 + (v(i,J,k) + v(i,J-1,k))**2) + grid_Ah = max(Ah(i,j), CS%min_grid_Ah) + grid_Re_Ah(i,j,k) = (sqrt(KE) * CS%grid_sp_h3(i,j)) / grid_Ah + enddo ; enddo + endif - endif ! biharmonic + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + d_del2u = G%IdyCu(I,j) * Del2u(I,j) - G%IdyCu(I-1,j) * Del2u(I-1,j) + d_del2v = G%IdxCv(i,J) * Del2v(i,J) - G%IdxCv(i,J-1) * Del2v(i,J-1) + d_str = Ah(i,j) * (CS%DY_dxT(i,j) * d_del2u - CS%DX_dyT(i,j) * d_del2v) - enddo ; enddo + str_xx(i,j) = str_xx(i,j) + d_str + + ! Keep a copy of the biharmonic contribution for backscatter parameterization + bhstr_xx(i,j) = d_str * (h(i,j,k) * CS%reduction_xx(i,j)) + enddo ; enddo + endif if (CS%biharmonic) then ! Gradient of Laplacian, for use in bi-harmonic term @@ -989,147 +1113,259 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, meke_res_fn = 1. + if ((CS%Smagorinsky_Kh) .or. (CS%Smagorinsky_Ah)) then + do J=js-1,Jeq ; do I=is-1,Ieq + sh_xy_sq = sh_xy(I,J) * sh_xy(I,J) + sh_xx_sq = 0.25 * ( & + (sh_xx(i,j) * sh_xx(i,j) + sh_xx(i+1,j+1) * sh_xx(i+1,j+1)) & + + (sh_xx(i,j+1) * sh_xx(i,j+1) + sh_xx(i+1,j) * sh_xx(i+1,j)) & + ) + Shear_mag(i,j) = sqrt(sh_xy_sq + sh_xx_sq) + enddo ; enddo + endif + do J=js-1,Jeq ; do I=is-1,Ieq - if ((CS%Smagorinsky_Kh) .or. (CS%Smagorinsky_Ah)) then - Shear_mag = sqrt(sh_xy(I,J)*sh_xy(I,J) + & - 0.25*((sh_xx(i,j)*sh_xx(i,j) + sh_xx(i+1,j+1)*sh_xx(i+1,j+1)) + & - (sh_xx(i,j+1)*sh_xx(i,j+1) + sh_xx(i+1,j)*sh_xx(i+1,j)))) + h2uq = 4.0 * (h_u(I,j) * h_u(I,j+1)) + h2vq = 4.0 * (h_v(i,J) * h_v(i+1,J)) + hq(I,J) = (2.0 * (h2uq * h2vq)) & + / (h_neglect3 + (h2uq + h2vq) * ((h_u(I,j) + h_u(I,j+1)) + (h_v(i,J) + h_v(i+1,J)))) + enddo ; enddo + + if (CS%better_bound_Ah .or. CS%better_bound_Kh) then + do J=js-1,Jeq ; do I=is-1,Ieq + h_min = min(h_u(I,j), h_u(I,j+1), h_v(i,J), h_v(i+1,J)) + hrat_min(i,j) = min(1.0, h_min / (hq(I,J) + h_neglect)) + enddo ; enddo + + if (CS%better_bound_Kh) then + do J=js-1,Jeq ; do I=is-1,Ieq + visc_bound_rem(i,j) = 1.0 + enddo ; enddo endif + endif + + if (CS%no_slip) then + do J=js-1,Jeq ; do I=is-1,Ieq + if (CS%no_slip .and. (G%mask2dBu(I,J) < 0.5)) then + if ((G%mask2dCu(I,j) + G%mask2dCu(I,j+1)) + & + (G%mask2dCv(i,J) + G%mask2dCv(i+1,J)) > 0.0) then + ! This is a coastal vorticity point, so modify hq and hrat_min. + + hu = G%mask2dCu(I,j) * h_u(I,j) + G%mask2dCu(I,j+1) * h_u(I,j+1) + hv = G%mask2dCv(i,J) * h_v(i,J) + G%mask2dCv(i+1,J) * h_v(i+1,J) + if ((G%mask2dCu(I,j) + G%mask2dCu(I,j+1)) * & + (G%mask2dCv(i,J) + G%mask2dCv(i+1,J)) == 0.0) then + ! Only one of hu and hv is nonzero, so just add them. + hq(I,J) = hu + hv + hrat_min(i,j) = 1.0 + else + ! Both hu and hv are nonzero, so take the harmonic mean. + hq(I,J) = 2.0 * (hu * hv) / ((hu + hv) + h_neglect) + hrat_min(i,j) = min(1.0, min(hu, hv) / (hq(I,J) + h_neglect) ) + endif + endif + endif + enddo ; enddo + endif + + if (CS%Laplacian) then if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then if (CS%use_QG_Leith_visc) then - vert_vort_mag = MIN(grad_vort_mag_q(I,J) + grad_div_mag_q(I,J), 3.*grad_vort_mag_q_2d(I,J)) + do J=js-1,Jeq ; do I=is-1,Ieq + grad_vort = grad_vort_mag_q(I,J) + grad_div_mag_q(I,J) + grad_vort_qg = 3. * grad_vort_mag_q_2d(I,J) + vert_vort_mag(i,j) = min(grad_vort, grad_vort_qg) + enddo ; enddo else - vert_vort_mag = (grad_vort_mag_q(I,J) + grad_div_mag_q(I,J)) + do J=js-1,Jeq ; do I=is-1,Ieq + vert_vort_mag(i,j) = grad_vort_mag_q(I,J) + grad_div_mag_q(I,J) + enddo ; enddo endif endif - h2uq = 4.0 * h_u(I,j) * h_u(I,j+1) - h2vq = 4.0 * h_v(i,J) * h_v(i+1,J) - hq(I,J) = 2.0 * h2uq * h2vq / (h_neglect3 + (h2uq + h2vq) * & - ((h_u(I,j) + h_u(I,j+1)) + (h_v(i,J) + h_v(i+1,J)))) - - if (CS%better_bound_Ah .or. CS%better_bound_Kh) then - hrat_min = min(1.0, min(h_u(I,j), h_u(I,j+1), h_v(i,J), h_v(i+1,J)) / & - (hq(I,J) + h_neglect) ) - visc_bound_rem = 1.0 - endif - if (CS%no_slip .and. (G%mask2dBu(I,J) < 0.5)) then - if ((G%mask2dCu(I,j) + G%mask2dCu(I,j+1)) + & - (G%mask2dCv(i,J) + G%mask2dCv(i+1,J)) > 0.0) then - ! This is a coastal vorticity point, so modify hq and hrat_min. - - hu = G%mask2dCu(I,j) * h_u(I,j) + G%mask2dCu(I,j+1) * h_u(I,j+1) - hv = G%mask2dCv(i,J) * h_v(i,J) + G%mask2dCv(i+1,J) * h_v(i+1,J) - if ((G%mask2dCu(I,j) + G%mask2dCu(I,j+1)) * & - (G%mask2dCv(i,J) + G%mask2dCv(i+1,J)) == 0.0) then - ! Only one of hu and hv is nonzero, so just add them. - hq(I,J) = hu + hv - hrat_min = 1.0 - else - ! Both hu and hv are nonzero, so take the harmonic mean. - hq(I,J) = 2.0 * (hu * hv) / ((hu + hv) + h_neglect) - hrat_min = min(1.0, min(hu, hv) / (hq(I,J) + h_neglect) ) - endif + ! Determine the Laplacian viscosity at q points, using the + ! largest value from several parameterizations. + + ! Static (pre-computed) background viscosity + do J=js-1,Jeq ; do I=is-1,Ieq + Kh(i,j) = CS%Kh_bg_xy(i,j) + enddo ; enddo + + if (CS%Smagorinsky_Kh) then + if (CS%add_LES_viscosity) then + do J=js-1,Jeq ; do I=is-1,Ieq + Kh(i,j) = Kh(i,j) + CS%Laplac2_const_xx(i,j) * Shear_mag(i,j) + enddo ; enddo + else + do J=js-1,Jeq ; do I=is-1,Ieq + Kh(i,j) = max(Kh(i,j), CS%Laplac2_const_xy(I,J) * Shear_mag(i,j) ) + enddo ; enddo endif endif - if (CS%Laplacian) then - ! Determine the Laplacian viscosity at q points, using the - ! largest value from several parameterizations. - Kh = CS%Kh_bg_xy(i,j) ! Static (pre-computed) background viscosity + if (CS%Leith_Kh) then if (CS%add_LES_viscosity) then - if (CS%Smagorinsky_Kh) Kh = Kh + CS%Laplac2_const_xx(i,j) * Shear_mag - if (CS%Leith_Kh) Kh = Kh + CS%Laplac3_const_xx(i,j) * vert_vort_mag*inv_PI3 + do J=js-1,Jeq ; do I=is-1,Ieq + Kh(i,j) = Kh(i,j) + CS%Laplac3_const_xx(i,j) * vert_vort_mag(i,j) * inv_PI3 + enddo ; enddo else - if (CS%Smagorinsky_Kh) Kh = max( Kh, CS%Laplac2_const_xy(I,J) * Shear_mag ) - if (CS%Leith_Kh) Kh = max( Kh, CS%Laplac3_const_xy(I,J) * vert_vort_mag*inv_PI3) + do J=js-1,Jeq ; do I=is-1,Ieq + Kh(i,j) = max(Kh(i,j), CS%Laplac3_const_xy(I,J) * vert_vort_mag(i,j) * inv_PI3) + enddo ; enddo endif - ! All viscosity contributions above are subject to resolution scaling - if (rescale_Kh) Kh = VarMix%Res_fn_q(i,j) * Kh - if (CS%res_scale_MEKE) meke_res_fn = VarMix%Res_fn_q(i,j) + endif + + ! All viscosity contributions above are subject to resolution scaling + + do J=js-1,Jeq ; do I=is-1,Ieq + ! NOTE: The following do-block can be decomposed and vectorized, but + ! appears to cause slowdown on some machines. Evidence suggests that + ! this is caused by excessive spilling of stack variables. + ! TODO: Vectorize these loops after stack usage has been reduced.. + + if (rescale_Kh) & + Kh(i,j) = VarMix%Res_fn_q(i,j) * Kh(i,j) + + if (CS%res_scale_MEKE) & + meke_res_fn = VarMix%Res_fn_q(i,j) + ! Older method of bounding for stability - if (legacy_bound) Kh = min(Kh, CS%Kh_Max_xy(i,j)) - Kh = max( Kh, CS%Kh_bg_min ) ! Place a floor on the viscosity, if desired. - if (use_MEKE_Ku) then ! *Add* the MEKE contribution (might be negative) - Kh = Kh + 0.25*( (MEKE%Ku(i,j) + MEKE%Ku(i+1,j+1)) + & + if (legacy_bound) & + Kh(i,j) = min(Kh(i,j), CS%Kh_Max_xy(i,j)) + + Kh(i,j) = max(Kh(i,j), CS%Kh_bg_min) ! Place a floor on the viscosity, if desired. + + if (use_MEKE_Ku) then + ! *Add* the MEKE contribution (might be negative) + Kh(i,j) = Kh(i,j) + 0.25*( (MEKE%Ku(i,j) + MEKE%Ku(i+1,j+1)) + & (MEKE%Ku(i+1,j) + MEKE%Ku(i,j+1)) ) * meke_res_fn endif + ! Older method of bounding for stability - if (CS%anisotropic) Kh = Kh + CS%Kh_aniso * CS%n1n2_q(I,J)**2 ! *Add* the shear component - ! of anisotropic viscosity + if (CS%anisotropic) & + ! *Add* the shear component of anisotropic viscosity + Kh(i,j) = Kh(i,j) + CS%Kh_aniso * CS%n1n2_q(I,J)**2 ! Newer method of bounding for stability if (CS%better_bound_Kh) then - if (Kh >= hrat_min*CS%Kh_Max_xy(I,J)) then - visc_bound_rem = 0.0 - Kh = hrat_min*CS%Kh_Max_xy(I,J) + if (Kh(i,j) >= hrat_min(i,j) * CS%Kh_Max_xy(I,J)) then + visc_bound_rem(i,j) = 0.0 + Kh(i,j) = hrat_min(i,j) * CS%Kh_Max_xy(I,J) elseif (CS%Kh_Max_xy(I,J)>0.) then - visc_bound_rem = 1.0 - Kh / (hrat_min*CS%Kh_Max_xy(I,J)) + visc_bound_rem(i,j) = 1.0 - Kh(i,j) / (hrat_min(i,j) * CS%Kh_Max_xy(I,J)) endif endif - if (CS%id_Kh_q>0 .or. CS%debug) Kh_q(I,J,k) = Kh - if (CS%id_vort_xy_q>0) vort_xy_q(I,J,k) = vort_xy(I,J) - if (CS%id_sh_xy_q>0) sh_xy_q(I,J,k) = sh_xy(I,J) + if (CS%id_Kh_q>0 .or. CS%debug) & + Kh_q(I,J,k) = Kh(i,j) - str_xy(I,J) = -Kh * sh_xy(I,J) - else ! not Laplacian - str_xy(I,J) = 0.0 - endif ! Laplacian + if (CS%id_vort_xy_q>0) & + vort_xy_q(I,J,k) = vort_xy(I,J) - if (CS%anisotropic) then + if (CS%id_sh_xy_q>0) & + sh_xy_q(I,J,k) = sh_xy(I,J) + enddo ; enddo + + do J=js-1,Jeq ; do I=is-1,Ieq + str_xy(I,J) = -Kh(i,j) * sh_xy(I,J) + enddo ; enddo + else + do J=js-1,Jeq ; do I=is-1,Ieq + str_xy(I,J) = 0. + enddo ; enddo + endif + + if (CS%anisotropic) then + do J=js-1,Jeq ; do I=is-1,Ieq ! Horizontal-tension averaged to q-points local_strain = 0.25 * ( (sh_xx(i,j) + sh_xx(i+1,j+1)) + (sh_xx(i+1,j) + sh_xx(i,j+1)) ) ! *Add* the tension contribution to the xy-component of stress str_xy(I,J) = str_xy(I,J) - CS%Kh_aniso * CS%n1n2_q(i,j) * CS%n1n1_m_n2n2_q(i,j) * local_strain - endif + enddo ; enddo + endif - if (CS%biharmonic) then + if (CS%biharmonic) then ! Determine the biharmonic viscosity at q points, using the ! largest value from several parameterizations. - AhSm = 0.0 ; AhLth = 0.0 - if (CS%Smagorinsky_Ah .or. CS%Leith_Ah) then - if (CS%Smagorinsky_Ah) then - if (CS%bound_Coriolis) then - AhSm = Shear_mag * (CS%Biharm_const_xy(I,J) + & - CS%Biharm_const2_xy(I,J)*Shear_mag) - else - AhSm = CS%Biharm_const_xy(I,J) * Shear_mag - endif - endif - if (CS%Leith_Ah) AhLth = CS%Biharm6_const_xy(I,J) * abs(Del2vort_q(I,J)) * inv_PI6 - Ah = MAX(MAX(CS%Ah_bg_xy(I,J), AhSm), AhLth) - if (CS%bound_Ah .and. .not.CS%better_bound_Ah) & - Ah = MIN(Ah, CS%Ah_Max_xy(I,J)) - else - Ah = CS%Ah_bg_xy(I,J) - endif ! Smagorinsky_Ah or Leith_Ah + do J=js-1,Jeq ; do I=is-1,Ieq + Ah(i,j) = CS%Ah_bg_xy(I,J) + enddo ; enddo - if (use_MEKE_Au) then ! *Add* the MEKE contribution - Ah = Ah + 0.25*( (MEKE%Au(i,j) + MEKE%Au(i+1,j+1)) + & - (MEKE%Au(i+1,j) + MEKE%Au(i,j+1)) ) + if (CS%Smagorinsky_Ah .or. CS%Leith_Ah) then + if (CS%Smagorinsky_Ah) then + if (CS%bound_Coriolis) then + do J=js-1,Jeq ; do I=is-1,Ieq + AhSm = Shear_mag(i,j) * (CS%Biharm_const_xy(I,J) & + + CS%Biharm_const2_xy(I,J) * Shear_mag(i,j) & + ) + Ah(i,j) = max(Ah(I,J), AhSm) + enddo ; enddo + else + do J=js-1,Jeq ; do I=is-1,Ieq + AhSm = CS%Biharm_const_xy(I,J) * Shear_mag(i,j) + Ah(i,j) = max(Ah(I,J), AhSm) + enddo ; enddo + endif endif - if (CS%Re_Ah > 0.0) then - KE = 0.125*((u(I,j,k)+u(I,j+1,k))**2 + (v(i,J,k)+v(i+1,J,k))**2) - Ah = sqrt(KE) * CS%Re_Ah_const_xy(i,j) + if (CS%Leith_Ah) then + do J=js-1,Jeq ; do I=is-1,Ieq + AhLth = CS%Biharm6_const_xy(I,J) * abs(Del2vort_q(I,J)) * inv_PI6 + Ah(i,j) = max(Ah(I,J), AhLth) + enddo ; enddo endif - if (CS%better_bound_Ah) then - Ah = MIN(Ah, visc_bound_rem*hrat_min*CS%Ah_Max_xy(I,J)) + if (CS%bound_Ah .and. .not.CS%better_bound_Ah) then + do J=js-1,Jeq ; do I=is-1,Ieq + Ah(i,j) = min(Ah(i,j), CS%Ah_Max_xy(I,J)) + enddo ; enddo endif + endif ! Smagorinsky_Ah or Leith_Ah - if (CS%id_Ah_q>0 .or. CS%debug) Ah_q(I,J,k) = Ah + if (use_MEKE_Au) then + ! *Add* the MEKE contribution + do J=js-1,Jeq ; do I=is-1,Ieq + Ah(i,j) = Ah(i,j) + 0.25 * ( & + (MEKE%Au(i,j) + MEKE%Au(i+1,j+1)) + (MEKE%Au(i+1,j) + MEKE%Au(i,j+1)) & + ) + enddo ; enddo + endif - str_xy(I,J) = str_xy(I,J) + Ah * ( dDel2vdx(I,J) + dDel2udy(I,J) ) + if (CS%Re_Ah > 0.0) then + do J=js-1,Jeq ; do I=is-1,Ieq + KE = 0.125 * ((u(I,j,k) + u(I,j+1,k))**2 + (v(i,J,k) + v(i+1,J,k))**2) + Ah(i,j) = sqrt(KE) * CS%Re_Ah_const_xy(i,j) + enddo ; enddo + endif - ! Keep a copy of the biharmonic contribution for backscatter parameterization - bhstr_xy(I,J) = Ah * ( dDel2vdx(I,J) + dDel2udy(I,J) ) * & - (hq(I,J) * G%mask2dBu(I,J) * CS%reduction_xy(I,J)) + if (CS%better_bound_Ah) then + if (CS%better_bound_Kh) then + do J=js-1,Jeq ; do I=is-1,Ieq + Ah(i,j) = min(Ah(i,j), visc_bound_rem(i,j) * hrat_min(i,j) * CS%Ah_Max_xy(I,J)) + enddo ; enddo + else + do J=js-1,Jeq ; do I=is-1,Ieq + Ah(i,j) = min(Ah(i,j), hrat_min(i,j) * CS%Ah_Max_xy(I,J)) + enddo ; enddo + endif + endif - endif ! biharmonic + if (CS%id_Ah_q>0 .or. CS%debug) then + do J=js-1,Jeq ; do I=is-1,Ieq + Ah_q(I,J,k) = Ah(i,j) + enddo ; enddo + endif - enddo ; enddo + ! Again, need to initialize str_xy as if its biharmonic + do J=js-1,Jeq ; do I=is-1,Ieq + d_str = Ah(i,j) * (dDel2vdx(I,J) + dDel2udy(I,J)) + + str_xy(I,J) = str_xy(I,J) + d_str + + ! Keep a copy of the biharmonic contribution for backscatter parameterization + bhstr_xy(I,J) = d_str * (hq(I,J) * G%mask2dBu(I,J) * CS%reduction_xy(I,J)) + enddo ; enddo + endif if (CS%use_GME) then call thickness_diffuse_get_KH(TD, KH_u_GME, KH_v_GME, G, GV) @@ -1197,14 +1433,15 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, str_xx(i,j) = str_xx(i,j) * (h(i,j,k) * CS%reduction_xx(i,j)) enddo ; enddo - do J=js-1,Jeq ; do I=is-1,Ieq - if (CS%no_slip) then + if (CS%no_slip) then + do J=js-1,Jeq ; do I=is-1,Ieq str_xy(I,J) = str_xy(I,J) * (hq(I,J) * CS%reduction_xy(I,J)) - else + enddo ; enddo + else + do J=js-1,Jeq ; do I=is-1,Ieq str_xy(I,J) = str_xy(I,J) * (hq(I,J) * G%mask2dBu(I,J) * CS%reduction_xy(I,J)) - endif - enddo ; enddo - + enddo ; enddo + endif endif ! use_GME ! Evaluate 1/h x.Div(h Grad u) or the biharmonic equivalent. @@ -1214,8 +1451,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, G%IdxCu(I,j)*(CS%dx2q(I,J-1)*str_xy(I,J-1) - & CS%dx2q(I,J) *str_xy(I,J))) * & G%IareaCu(I,j)) / (h_u(I,j) + h_neglect) - enddo ; enddo + if (apply_OBC) then ! This is not the right boundary condition. If all the masking of tendencies are done ! correctly later then eliminating this block should not change answers. @@ -1237,6 +1474,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, CS%dx2h(i,j+1)*str_xx(i,j+1))) * & G%IareaCv(i,J)) / (h_v(i,J) + h_neglect) enddo ; enddo + if (apply_OBC) then ! This is not the right boundary condition. If all the masking of tendencies are done ! correctly later then eliminating this block should not change answers. @@ -1288,28 +1526,27 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, do j=js,je ; do i=is,ie FatH = 0.25*( (abs(G%CoriolisBu(I-1,J-1)) + abs(G%CoriolisBu(I,J))) + & (abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J-1))) ) - Shear_mag = sqrt(sh_xx(i,j)*sh_xx(i,j) + & + Shear_mag_bc = sqrt(sh_xx(i,j) * sh_xx(i,j) + & 0.25*((sh_xy(I-1,J-1)*sh_xy(I-1,J-1) + sh_xy(I,J)*sh_xy(I,J)) + & (sh_xy(I-1,J)*sh_xy(I-1,J) + sh_xy(I,J-1)*sh_xy(I,J-1)))) if (CS%answers_2018) then FatH = (US%s_to_T*FatH)**MEKE%backscatter_Ro_pow ! f^n ! Note the hard-coded dimensional constant in the following line that can not ! be rescaled for dimensional consistency. - Shear_mag = ( ( (US%s_to_T*Shear_mag)**MEKE%backscatter_Ro_pow ) + 1.e-30 ) & + Shear_mag_bc = (((US%s_to_T * Shear_mag_bc)**MEKE%backscatter_Ro_pow) + 1.e-30) & * MEKE%backscatter_Ro_c ! c * D^n ! The Rossby number function is g(Ro) = 1/(1+c.Ro^n) ! RoScl = 1 - g(Ro) - RoScl = Shear_mag / ( FatH + Shear_mag ) ! = 1 - f^n/(f^n+c*D^n) + RoScl = Shear_mag_bc / (FatH + Shear_mag_bc) ! = 1 - f^n/(f^n+c*D^n) else - if (FatH <= backscat_subround*Shear_mag) then + if (FatH <= backscat_subround*Shear_mag_bc) then RoScl = 1.0 else - Sh_F_pow = MEKE%backscatter_Ro_c * (Shear_mag / FatH)**MEKE%backscatter_Ro_pow + Sh_F_pow = MEKE%backscatter_Ro_c * (Shear_mag_bc / FatH)**MEKE%backscatter_Ro_pow RoScl = Sh_F_pow / (1.0 + Sh_F_pow) ! = 1 - f^n/(f^n+c*D^n) endif endif - MEKE%mom_src(i,j) = MEKE%mom_src(i,j) + GV%H_to_RZ * ( & ((str_xx(i,j)-RoScl*bhstr_xx(i,j))*(u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j) & -(str_xx(i,j)-RoScl*bhstr_xx(i,j))*(v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j)) & From 571013dad0a7f4971f629cbda16b1c396620d812 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 12 Jan 2021 15:32:40 -0500 Subject: [PATCH 03/11] +Added MOM_coms_wrapper.F90 Added the new module MOM_coms_wrapper, along with explicit interfaces for the broadcast routine for the cases that might actually be used by MOM6. With these new interfaces, the source PE has been made into an optional argument, and there is a new optional argument to indicate whether the broadcast is blocking. Also the MOM_horizontal_regridding module has been updated to reflect these changes. All answers are bitwise identical, but an existing required argument to broadcast has been made optional and there is a new optional argument. --- src/framework/MOM_coms.F90 | 19 +-- src/framework/MOM_coms_wrapper.F90 | 160 ++++++++++++++++++++ src/framework/MOM_horizontal_regridding.F90 | 52 +++---- 3 files changed, 183 insertions(+), 48 deletions(-) create mode 100644 src/framework/MOM_coms_wrapper.F90 diff --git a/src/framework/MOM_coms.F90 b/src/framework/MOM_coms.F90 index 04ed46ad22..13fc4df75d 100644 --- a/src/framework/MOM_coms.F90 +++ b/src/framework/MOM_coms.F90 @@ -5,22 +5,19 @@ module MOM_coms ! This file is part of MOM6. See LICENSE.md for the license. use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING -use fms_mod, only : fms_end, MOM_infra_init => fms_init -use memutils_mod, only : print_memuse_stats -use mpp_mod, only : PE_here => mpp_pe, root_PE => mpp_root_pe, num_PEs => mpp_npes -use mpp_mod, only : Set_PElist => mpp_set_current_pelist, Get_PElist => mpp_get_current_pelist -use mpp_mod, only : broadcast => mpp_broadcast, field_chksum => mpp_chksum -use mpp_mod, only : sum_across_PEs => mpp_sum, max_across_PEs => mpp_max, min_across_PEs => mpp_min +use MOM_coms_wrapper, only : PE_here, root_PE, num_PEs, Set_PElist, Get_PElist +use MOM_coms_wrapper, only : broadcast, field_chksum, MOM_infra_init, MOM_infra_end +use MOM_coms_wrapper, only : sum_across_PEs, max_across_PEs, min_across_PEs implicit none ; private public :: PE_here, root_PE, num_PEs, MOM_infra_init, MOM_infra_end public :: broadcast, sum_across_PEs, min_across_PEs, max_across_PEs, field_chksum +public :: Set_PElist, Get_PElist public :: reproducing_sum, reproducing_sum_EFP, EFP_sum_across_PEs, EFP_list_sum_across_PEs public :: EFP_plus, EFP_minus, EFP_to_real, real_to_EFP, EFP_real_diff public :: operator(+), operator(-), assignment(=) public :: query_EFP_overflow_error, reset_EFP_overflow_error -public :: Set_PElist, Get_PElist ! This module provides interfaces to the non-domain-oriented communication subroutines. @@ -880,12 +877,4 @@ subroutine EFP_val_sum_across_PEs(EFP, error) end subroutine EFP_val_sum_across_PEs - -!> This subroutine carries out all of the calls required to close out the infrastructure cleanly. -!! This should only be called in ocean-only runs, as the coupler takes care of this in coupled runs. -subroutine MOM_infra_end - call print_memuse_stats( 'Memory HiWaterMark', always=.TRUE. ) - call fms_end -end subroutine MOM_infra_end - end module MOM_coms diff --git a/src/framework/MOM_coms_wrapper.F90 b/src/framework/MOM_coms_wrapper.F90 new file mode 100644 index 0000000000..954f6da93c --- /dev/null +++ b/src/framework/MOM_coms_wrapper.F90 @@ -0,0 +1,160 @@ +!> Thin interfaces to non-domain-oriented mpp communication subroutines +module MOM_coms_wrapper + +! This file is part of MOM6. See LICENSE.md for the license. + +use fms_mod, only : fms_end, MOM_infra_init => fms_init +use memutils_mod, only : print_memuse_stats +use mpp_mod, only : PE_here => mpp_pe, root_PE => mpp_root_pe, num_PEs => mpp_npes +use mpp_mod, only : Set_PElist => mpp_set_current_pelist, Get_PElist => mpp_get_current_pelist +use mpp_mod, only : mpp_broadcast, mpp_sync, mpp_sync_self, field_chksum => mpp_chksum +use mpp_mod, only : sum_across_PEs => mpp_sum, max_across_PEs => mpp_max, min_across_PEs => mpp_min + +implicit none ; private + +public :: PE_here, root_PE, num_PEs, MOM_infra_init, MOM_infra_end, Set_PElist, Get_PElist +public :: broadcast, sum_across_PEs, min_across_PEs, max_across_PEs, field_chksum + +! This module provides interfaces to the non-domain-oriented communication subroutines. + +!> Communicate an array, string or scalar from one PE to others +interface broadcast + module procedure broadcast_char, broadcast_int0D, broadcast_int1D + module procedure broadcast_real0D, broadcast_real1D, broadcast_real2D +end interface broadcast + +contains + +!> Communicate a 1-D array of character strings from one PE to others +subroutine broadcast_char(dat, length, from_PE, PElist, blocking) + character(len=*), intent(inout) :: dat(:) !< The data to communicate and destination + integer, intent(in) :: length !< The length of each string + integer, optional, intent(in) :: from_PE !< The source PE, by default the root PE + integer, optional, intent(in) :: PElist(:) !< The list of participating PEs, by default the + !! active PE set as previously set via Set_PElist. + logical, optional, intent(in) :: blocking !< If true, barriers are added around the call + + integer :: src_PE ! The processor that is sending the data + logical :: do_block ! If true add synchronizing barriers + + do_block = .false. ; if (present(blocking)) do_block = blocking + if (present(from_PE)) then ; src_PE = from_PE ; else ; src_PE = root_PE() ; endif + + if (do_block) call mpp_sync(PElist) + call mpp_broadcast(dat, length, src_PE, PElist) + if (do_block) call mpp_sync_self(PElist) + +end subroutine broadcast_char + +!> Communicate an integer from one PE to others +subroutine broadcast_int0D(dat, from_PE, PElist, blocking) + integer, intent(inout) :: dat !< The data to communicate and destination + integer, optional, intent(in) :: from_PE !< The source PE, by default the root PE + integer, optional, intent(in) :: PElist(:) !< The list of participating PEs, by default the + !! active PE set as previously set via Set_PElist. + logical, optional, intent(in) :: blocking !< If true, barriers are added around the call + + integer :: src_PE ! The processor that is sending the data + logical :: do_block ! If true add synchronizing barriers + + do_block = .false. ; if (present(blocking)) do_block = blocking + if (present(from_PE)) then ; src_PE = from_PE ; else ; src_PE = root_PE() ; endif + + if (do_block) call mpp_sync(PElist) + call mpp_broadcast(dat, src_PE, PElist) + if (do_block) call mpp_sync_self(PElist) + +end subroutine broadcast_int0D + +!> Communicate a 1-D array of integers from one PE to others +subroutine broadcast_int1D(dat, length, from_PE, PElist, blocking) + integer, dimension(:), intent(inout) :: dat !< The data to communicate and destination + integer, intent(in) :: length !< The number of data elements + integer, optional, intent(in) :: from_PE !< The source PE, by default the root PE + integer, optional, intent(in) :: PElist(:) !< The list of participating PEs, by default the + !! active PE set as previously set via Set_PElist. + logical, optional, intent(in) :: blocking !< If true, barriers are added around the call + + integer :: src_PE ! The processor that is sending the data + logical :: do_block ! If true add synchronizing barriers + + do_block = .false. ; if (present(blocking)) do_block = blocking + if (present(from_PE)) then ; src_PE = from_PE ; else ; src_PE = root_PE() ; endif + + if (do_block) call mpp_sync(PElist) + call mpp_broadcast(dat, length, src_PE, PElist) + if (do_block) call mpp_sync_self(PElist) + +end subroutine broadcast_int1D + +!> Communicate a real number from one PE to others +subroutine broadcast_real0D(dat, from_PE, PElist, blocking) + real, intent(inout) :: dat !< The data to communicate and destination + integer, optional, intent(in) :: from_PE !< The source PE, by default the root PE + integer, optional, intent(in) :: PElist(:) !< The list of participating PEs, by default the + !! active PE set as previously set via Set_PElist. + logical, optional, intent(in) :: blocking !< If true, barriers are added around the call + + integer :: src_PE ! The processor that is sending the data + logical :: do_block ! If true add synchronizing barriers + + do_block = .false. ; if (present(blocking)) do_block = blocking + if (present(from_PE)) then ; src_PE = from_PE ; else ; src_PE = root_PE() ; endif + + if (do_block) call mpp_sync(PElist) + call mpp_broadcast(dat, src_PE, PElist) + if (do_block) call mpp_sync_self(PElist) + +end subroutine broadcast_real0D + +!> Communicate a 1-D array of reals from one PE to others +subroutine broadcast_real1D(dat, length, from_PE, PElist, blocking) + real, dimension(:), intent(inout) :: dat !< The data to communicate and destination + integer, intent(in) :: length !< The number of data elements + integer, optional, intent(in) :: from_PE !< The source PE, by default the root PE + integer, optional, intent(in) :: PElist(:) !< The list of participating PEs, by default the + !! active PE set as previously set via Set_PElist. + logical, optional, intent(in) :: blocking !< If true, barriers are added around the call + + integer :: src_PE ! The processor that is sending the data + logical :: do_block ! If true add synchronizing barriers + + do_block = .false. ; if (present(blocking)) do_block = blocking + if (present(from_PE)) then ; src_PE = from_PE ; else ; src_PE = root_PE() ; endif + + if (do_block) call mpp_sync(PElist) + call mpp_broadcast(dat, length, src_PE, PElist) + if (do_block) call mpp_sync_self(PElist) + +end subroutine broadcast_real1D + +!> Communicate a 2-D array of reals from one PE to others +subroutine broadcast_real2D(dat, length, from_PE, PElist, blocking) + real, dimension(:,:), intent(inout) :: dat !< The data to communicate and destination + integer, intent(in) :: length !< The total number of data elements + integer, optional, intent(in) :: from_PE !< The source PE, by default the root PE + integer, optional, intent(in) :: PElist(:) !< The list of participating PEs, by default the + !! active PE set as previously set via Set_PElist. + logical, optional, intent(in) :: blocking !< If true, barriers are added around the call + + integer :: src_PE ! The processor that is sending the data + logical :: do_block ! If true add synchronizing barriers + + do_block = .false. ; if (present(blocking)) do_block = blocking + if (present(from_PE)) then ; src_PE = from_PE ; else ; src_PE = root_PE() ; endif + + if (do_block) call mpp_sync(PElist) + call mpp_broadcast(dat, length, src_PE, PElist) + if (do_block) call mpp_sync_self(PElist) + +end subroutine broadcast_real2D + + +!> This subroutine carries out all of the calls required to close out the infrastructure cleanly. +!! This should only be called in ocean-only runs, as the coupler takes care of this in coupled runs. +subroutine MOM_infra_end + call print_memuse_stats( 'Memory HiWaterMark', always=.TRUE. ) + call fms_end() +end subroutine MOM_infra_end + +end module MOM_coms_wrapper diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index 4f98038f12..9b340f3aa7 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -3,32 +3,22 @@ module MOM_horizontal_regridding ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_debugging, only : hchksum -use MOM_coms, only : max_across_PEs, min_across_PEs -use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end -use MOM_cpu_clock, only : CLOCK_ROUTINE, CLOCK_LOOP -use MOM_domains, only : pass_var, pass_vector, sum_across_PEs, broadcast -use MOM_domains, only : root_PE, To_All, SCALAR_PAIR, CGRID_NE, AGRID +use MOM_debugging, only : hchksum +use MOM_coms, only : max_across_PEs, min_across_PEs, sum_across_PEs, broadcast +use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_LOOP +use MOM_domains, only : pass_var use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint -use MOM_file_parser, only : get_param, read_param, log_param, param_file_type -use MOM_file_parser, only : log_version -use MOM_get_input, only : directories -use MOM_grid, only : ocean_grid_type, isPointInCell -use MOM_io, only : close_file, fieldtype, file_exists -use MOM_io, only : open_file, read_data, read_axis_data, SINGLE_FILE, MULTIPLE -use MOM_io, only : slasher, vardesc, write_field -use MOM_string_functions, only : uppercase -use MOM_time_manager, only : time_type, get_external_field_size -use MOM_time_manager, only : init_external_field -use MOM_time_manager, only : get_external_field_axes, get_external_field_missing +use MOM_file_parser, only : get_param, log_param, log_version, param_file_type +use MOM_grid, only : ocean_grid_type +use MOM_io_wrapper, only : axistype, get_axis_data +use MOM_time_manager, only : time_type +use MOM_time_manager, only : init_external_field, get_external_field_size +use MOM_time_manager, only : get_external_field_axes, get_external_field_missing use MOM_transform_FMS, only : time_interp_external => rotated_time_interp_external -use MOM_variables, only : thermo_var_ptrs -use mpp_io_mod, only : axistype, mpp_get_axis_data -use mpp_mod, only : mpp_broadcast, mpp_sync, mpp_sync_self, mpp_max -use horiz_interp_mod, only : horiz_interp_new, horiz_interp, horiz_interp_type -use horiz_interp_mod, only : horiz_interp_init, horiz_interp_del +use horiz_interp_mod, only : horiz_interp_new, horiz_interp, horiz_interp_type +use horiz_interp_mod, only : horiz_interp_init, horiz_interp_del use netcdf @@ -463,7 +453,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, endif max_depth = maxval(G%bathyT) - call mpp_max(max_depth) + call max_across_PEs(max_depth) if (z_edges_in(kd+1) Date: Wed, 13 Jan 2021 07:10:45 -0500 Subject: [PATCH 04/11] +Added MOM_interpolate.F90 Added the new module MOM_interpolate to wrap the time_interp_external and horiz_interp modules. The overloaded interface time_interp_external had to be renamed to time_interp_extern in the version offered by MOM_interpolate because of a weird compile time problem with the PGI 19.10.0 compiler. Some large blocks of unused code in MOM_horizontal_regridding were commented out. Also modified 6 files outside of the framework directory to use these new interfaces, but they will all work without these changes. All answers are bitwise identical. --- .../MOM_surface_forcing_gfdl.F90 | 49 ++-- config_src/coupled_driver/ocean_model_MOM.F90 | 22 +- config_src/solo_driver/MOM_driver.F90 | 14 +- src/core/MOM_open_boundary.F90 | 9 +- src/framework/MOM_horizontal_regridding.F90 | 234 +++++++++--------- src/framework/MOM_interpolate.F90 | 143 +++++++++++ src/ice_shelf/MOM_ice_shelf.F90 | 14 +- .../vertical/MOM_diabatic_aux.F90 | 8 +- 8 files changed, 313 insertions(+), 180 deletions(-) create mode 100644 src/framework/MOM_interpolate.F90 diff --git a/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 b/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 index be960accd6..67f6643a42 100644 --- a/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 +++ b/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 @@ -5,7 +5,7 @@ module MOM_surface_forcing_gfdl !#CTRL# use MOM_controlled_forcing, only : apply_ctrl_forcing, register_ctrl_forcing_restarts !#CTRL# use MOM_controlled_forcing, only : controlled_forcing_init, controlled_forcing_end !#CTRL# use MOM_controlled_forcing, only : ctrl_forcing_CS -use MOM_coms, only : reproducing_sum +use MOM_coms, only : reproducing_sum, field_chksum use MOM_constants, only : hlv, hlf use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_SUBCOMPONENT @@ -21,6 +21,8 @@ module MOM_surface_forcing_gfdl use MOM_forcing_type, only : allocate_mech_forcing, deallocate_mech_forcing use MOM_get_input, only : Get_MOM_Input, directories use MOM_grid, only : ocean_grid_type +use MOM_interpolate, only : init_external_field, time_interp_extern +use MOM_interpolate, only : time_interp_external_init use MOM_io, only : slasher, write_version_number, MOM_read_data use MOM_restart, only : register_restart_field, restart_init, MOM_restart_CS use MOM_restart, only : restart_init_end, save_restart, restore_state @@ -36,9 +38,6 @@ module MOM_surface_forcing_gfdl use coupler_types_mod, only : coupler_type_copy_data use data_override_mod, only : data_override_init, data_override use fms_mod, only : stdout -use mpp_mod, only : mpp_chksum -use time_interp_external_mod, only : init_external_field, time_interp_external -use time_interp_external_mod, only : time_interp_external_init implicit none ; private @@ -350,7 +349,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! Salinity restoring logic if (CS%restore_salt) then - call time_interp_external(CS%id_srestore,Time,data_restore) + call time_interp_extern(CS%id_srestore, Time, data_restore) ! open_ocn_mask indicates where to restore salinity (1 means restore, 0 does not) open_ocn_mask(:,:) = 1.0 if (CS%mask_srestore_under_ice) then ! Do not restore under sea-ice @@ -407,7 +406,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! SST restoring logic if (CS%restore_temp) then - call time_interp_external(CS%id_trestore,Time,data_restore) + call time_interp_extern(CS%id_trestore, Time, data_restore) do j=js,je ; do i=is,ie delta_sst = data_restore(i,j)- sfc_state%SST(i,j) delta_sst = sign(1.0,delta_sst)*min(abs(delta_sst),CS%max_delta_trestore) @@ -1486,7 +1485,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, wind_stagger) enddo ; enddo endif - call time_interp_external_init + call time_interp_external_init() ! Optionally read a x-y gustiness field in place of a global constant. call get_param(param_file, mdl, "READ_GUST_2D", CS%read_gust_2d, & @@ -1632,27 +1631,27 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt) outunit = stdout() write(outunit,*) "BEGIN CHECKSUM(ice_ocean_boundary_type):: ", id, timestep - write(outunit,100) 'iobt%u_flux ', mpp_chksum( iobt%u_flux ) - write(outunit,100) 'iobt%v_flux ', mpp_chksum( iobt%v_flux ) - write(outunit,100) 'iobt%t_flux ', mpp_chksum( iobt%t_flux ) - write(outunit,100) 'iobt%q_flux ', mpp_chksum( iobt%q_flux ) - write(outunit,100) 'iobt%salt_flux ', mpp_chksum( iobt%salt_flux ) - write(outunit,100) 'iobt%lw_flux ', mpp_chksum( iobt%lw_flux ) - write(outunit,100) 'iobt%sw_flux_vis_dir', mpp_chksum( iobt%sw_flux_vis_dir) - write(outunit,100) 'iobt%sw_flux_vis_dif', mpp_chksum( iobt%sw_flux_vis_dif) - write(outunit,100) 'iobt%sw_flux_nir_dir', mpp_chksum( iobt%sw_flux_nir_dir) - write(outunit,100) 'iobt%sw_flux_nir_dif', mpp_chksum( iobt%sw_flux_nir_dif) - write(outunit,100) 'iobt%lprec ', mpp_chksum( iobt%lprec ) - write(outunit,100) 'iobt%fprec ', mpp_chksum( iobt%fprec ) - write(outunit,100) 'iobt%runoff ', mpp_chksum( iobt%runoff ) - write(outunit,100) 'iobt%calving ', mpp_chksum( iobt%calving ) - write(outunit,100) 'iobt%p ', mpp_chksum( iobt%p ) + write(outunit,100) 'iobt%u_flux ', field_chksum( iobt%u_flux ) + write(outunit,100) 'iobt%v_flux ', field_chksum( iobt%v_flux ) + write(outunit,100) 'iobt%t_flux ', field_chksum( iobt%t_flux ) + write(outunit,100) 'iobt%q_flux ', field_chksum( iobt%q_flux ) + write(outunit,100) 'iobt%salt_flux ', field_chksum( iobt%salt_flux ) + write(outunit,100) 'iobt%lw_flux ', field_chksum( iobt%lw_flux ) + write(outunit,100) 'iobt%sw_flux_vis_dir', field_chksum( iobt%sw_flux_vis_dir) + write(outunit,100) 'iobt%sw_flux_vis_dif', field_chksum( iobt%sw_flux_vis_dif) + write(outunit,100) 'iobt%sw_flux_nir_dir', field_chksum( iobt%sw_flux_nir_dir) + write(outunit,100) 'iobt%sw_flux_nir_dif', field_chksum( iobt%sw_flux_nir_dif) + write(outunit,100) 'iobt%lprec ', field_chksum( iobt%lprec ) + write(outunit,100) 'iobt%fprec ', field_chksum( iobt%fprec ) + write(outunit,100) 'iobt%runoff ', field_chksum( iobt%runoff ) + write(outunit,100) 'iobt%calving ', field_chksum( iobt%calving ) + write(outunit,100) 'iobt%p ', field_chksum( iobt%p ) if (associated(iobt%ustar_berg)) & - write(outunit,100) 'iobt%ustar_berg ', mpp_chksum( iobt%ustar_berg ) + write(outunit,100) 'iobt%ustar_berg ', field_chksum( iobt%ustar_berg ) if (associated(iobt%area_berg)) & - write(outunit,100) 'iobt%area_berg ', mpp_chksum( iobt%area_berg ) + write(outunit,100) 'iobt%area_berg ', field_chksum( iobt%area_berg ) if (associated(iobt%mass_berg)) & - write(outunit,100) 'iobt%mass_berg ', mpp_chksum( iobt%mass_berg ) + write(outunit,100) 'iobt%mass_berg ', field_chksum( iobt%mass_berg ) 100 FORMAT(" CHECKSUM::",A20," = ",Z20) call coupler_type_write_chksums(iobt%fluxes, outunit, 'iobt%') diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index b429da649b..12f803a970 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -15,6 +15,7 @@ module ocean_model_mod use MOM, only : extract_surface_state, allocate_surface_state, finish_MOM_initialization use MOM, only : get_MOM_state_elements, MOM_state_is_synchronized use MOM, only : get_ocean_stocks, step_offline +use MOM_coms, only : field_chksum use MOM_constants, only : CELSIUS_KELVIN_OFFSET, hlf use MOM_diag_mediator, only : diag_ctrl, enable_averaging, disable_averaging use MOM_diag_mediator, only : diag_mediator_close_registration, diag_mediator_end @@ -22,6 +23,7 @@ module ocean_model_mod use MOM_domains, only : TO_ALL, Omit_Corners use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave +use MOM_EOS, only : gsw_sp_from_sr, gsw_pt_from_ct use MOM_file_parser, only : get_param, log_version, close_param_file, param_file_type use MOM_forcing_type, only : forcing, mech_forcing, allocate_forcing_type use MOM_forcing_type, only : fluxes_accumulate, get_net_mass_forcing @@ -48,6 +50,8 @@ module ocean_model_mod use MOM_verticalGrid, only : verticalGrid_type use MOM_ice_shelf, only : initialize_ice_shelf, shelf_calc_flux, ice_shelf_CS use MOM_ice_shelf, only : add_shelf_forces, ice_shelf_end, ice_shelf_save_restart +use MOM_wave_interface, only: wave_parameters_CS, MOM_wave_interface_init +use MOM_wave_interface, only: MOM_wave_interface_init_lite, Update_Surface_Waves use coupler_types_mod, only : coupler_1d_bc_type, coupler_2d_bc_type use coupler_types_mod, only : coupler_type_spawn, coupler_type_write_chksums use coupler_types_mod, only : coupler_type_initialized, coupler_type_copy_data @@ -56,10 +60,6 @@ module ocean_model_mod use mpp_domains_mod, only : mpp_define_domains, mpp_get_compute_domain, mpp_get_data_domain use atmos_ocean_fluxes_mod, only : aof_set_coupler_flux use fms_mod, only : stdout -use mpp_mod, only : mpp_chksum -use MOM_EOS, only : gsw_sp_from_sr, gsw_pt_from_ct -use MOM_wave_interface, only: wave_parameters_CS, MOM_wave_interface_init -use MOM_wave_interface, only: MOM_wave_interface_init_lite, Update_Surface_Waves #include @@ -1130,13 +1130,13 @@ subroutine ocean_public_type_chksum(id, timestep, ocn) outunit = stdout() write(outunit,*) "BEGIN CHECKSUM(ocean_type):: ", id, timestep - write(outunit,100) 'ocean%t_surf ',mpp_chksum(ocn%t_surf ) - write(outunit,100) 'ocean%s_surf ',mpp_chksum(ocn%s_surf ) - write(outunit,100) 'ocean%u_surf ',mpp_chksum(ocn%u_surf ) - write(outunit,100) 'ocean%v_surf ',mpp_chksum(ocn%v_surf ) - write(outunit,100) 'ocean%sea_lev ',mpp_chksum(ocn%sea_lev) - write(outunit,100) 'ocean%frazil ',mpp_chksum(ocn%frazil ) - write(outunit,100) 'ocean%melt_potential ',mpp_chksum(ocn%melt_potential) + write(outunit,100) 'ocean%t_surf ', field_chksum(ocn%t_surf ) + write(outunit,100) 'ocean%s_surf ', field_chksum(ocn%s_surf ) + write(outunit,100) 'ocean%u_surf ', field_chksum(ocn%u_surf ) + write(outunit,100) 'ocean%v_surf ', field_chksum(ocn%v_surf ) + write(outunit,100) 'ocean%sea_lev ', field_chksum(ocn%sea_lev) + write(outunit,100) 'ocean%frazil ', field_chksum(ocn%frazil ) + write(outunit,100) 'ocean%melt_potential ', field_chksum(ocn%melt_potential) call coupler_type_write_chksums(ocn%fields, outunit, 'ocean%') 100 FORMAT(" CHECKSUM::",A20," = ",Z20) diff --git a/config_src/solo_driver/MOM_driver.F90 b/config_src/solo_driver/MOM_driver.F90 index 9726aa1281..c9383a4287 100644 --- a/config_src/solo_driver/MOM_driver.F90 +++ b/config_src/solo_driver/MOM_driver.F90 @@ -32,6 +32,7 @@ program MOM_main use MOM, only : extract_surface_state, finish_MOM_initialization use MOM, only : get_MOM_state_elements, MOM_state_is_synchronized use MOM, only : step_offline + use MOM_coms, only : Set_PElist use MOM_domains, only : MOM_infra_init, MOM_infra_end use MOM_error_handler, only : MOM_error, MOM_mesg, WARNING, FATAL, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint @@ -41,6 +42,7 @@ program MOM_main use MOM_forcing_type, only : mech_forcing_diags, MOM_forcing_chksum, MOM_mech_forcing_chksum use MOM_get_input, only : directories use MOM_grid, only : ocean_grid_type + use MOM_interpolate, only : time_interp_external_init use MOM_io, only : file_exists, open_file, close_file use MOM_io, only : check_nml_error, io_infra_init, io_infra_end use MOM_io, only : APPEND_FILE, ASCII_FILE, READONLY_FILE, SINGLE_FILE @@ -64,8 +66,6 @@ program MOM_main use MOM_get_input, only : get_MOM_input use ensemble_manager_mod, only : ensemble_manager_init, get_ensemble_size use ensemble_manager_mod, only : ensemble_pelist_setup - use mpp_mod, only : set_current_pelist => mpp_set_current_pelist - use time_interp_external_mod, only : time_interp_external_init use fms_affinity_mod, only : fms_affinity_init, fms_affinity_set,fms_affinity_get use MOM_ice_shelf, only : initialize_ice_shelf, ice_shelf_end, ice_shelf_CS @@ -229,7 +229,7 @@ program MOM_main allocate(ocean_pelist(nPEs_per)) call ensemble_pelist_setup(.true., 0, nPEs_per, 0, 0, atm_pelist, ocean_pelist, & land_pelist, ice_pelist) - call set_current_pelist(ocean_pelist) + call Set_PElist(ocean_pelist) deallocate(ocean_pelist) endif @@ -286,17 +286,17 @@ program MOM_main if (sum(date_init) > 0) then - Start_time = set_date(date_init(1),date_init(2), date_init(3), & - date_init(4),date_init(5),date_init(6)) + Start_time = set_date(date_init(1), date_init(2), date_init(3), & + date_init(4), date_init(5), date_init(6)) else Start_time = real_to_time(0.0) endif - call time_interp_external_init + call time_interp_external_init() if (sum(date) >= 0) then ! In this case, the segment starts at a time fixed by ocean_solo.res - segment_start_time = set_date(date(1),date(2),date(3),date(4),date(5),date(6)) + segment_start_time = set_date(date(1), date(2), date(3), date(4), date(5), date(6)) Time = segment_start_time else ! In this case, the segment starts at a time read from the MOM restart file diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 3fc0c9bcba..4faab21177 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -24,8 +24,7 @@ module MOM_open_boundary use MOM_tidal_forcing, only : astro_longitudes, astro_longitudes_init, eq_phase, nodal_fu, tidal_frequency use MOM_time_manager, only : set_date, time_type, time_type_to_real, operator(-) use MOM_tracer_registry, only : tracer_type, tracer_registry_type, tracer_name_lookup -use time_interp_external_mod, only : init_external_field, time_interp_external -use time_interp_external_mod, only : time_interp_external_init +use MOM_interpolate, only : init_external_field, time_interp_extern, time_interp_external_init use MOM_remapping, only : remappingSchemesDoc, remappingDefaultScheme, remapping_CS use MOM_remapping, only : initialize_remapping, remapping_core_h, end_remapping use MOM_regridding, only : regridding_CS @@ -3897,7 +3896,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) ! TODO: Since we conditionally rotate a subset of tmp_buffer_in after ! reading the value, it is currently not possible to use the rotated - ! implementation of time_interp_external. + ! implementation of time_interp_extern. ! For now, we must explicitly allocate and rotate this array. if (turns /= 0) then if (modulo(turns, 2) /= 0) then @@ -3909,7 +3908,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) tmp_buffer_in => tmp_buffer endif - call time_interp_external(segment%field(m)%fid,Time, tmp_buffer_in) + call time_interp_extern(segment%field(m)%fid,Time, tmp_buffer_in) ! NOTE: Rotation of face-points require that we skip the final value if (turns /= 0) then ! TODO: This is hardcoded for 90 degrees, and needs to be generalized. @@ -3976,7 +3975,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) ! no dz for tidal variables if (segment%field(m)%nk_src > 1 .and.& (index(segment%field(m)%name, 'phase') .le. 0 .and. index(segment%field(m)%name, 'amp') .le. 0)) then - call time_interp_external(segment%field(m)%fid_dz,Time, tmp_buffer_in) + call time_interp_extern(segment%field(m)%fid_dz,Time, tmp_buffer_in) if (turns /= 0) then ! TODO: This is hardcoded for 90 degrees, and needs to be generalized. if (segment%is_E_or_W & diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index 9b340f3aa7..6e09d391d9 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -13,12 +13,8 @@ module MOM_horizontal_regridding use MOM_grid, only : ocean_grid_type use MOM_io_wrapper, only : axistype, get_axis_data use MOM_time_manager, only : time_type -use MOM_time_manager, only : init_external_field, get_external_field_size -use MOM_time_manager, only : get_external_field_axes, get_external_field_missing -use MOM_transform_FMS, only : time_interp_external => rotated_time_interp_external - -use horiz_interp_mod, only : horiz_interp_new, horiz_interp, horiz_interp_type -use horiz_interp_mod, only : horiz_interp_init, horiz_interp_del +use MOM_interpolate, only : time_interp_extern, get_external_field_info, horiz_interp_init +use MOM_interpolate, only : horiz_interp_new, horiz_interp, horiz_interp_type use netcdf @@ -31,10 +27,10 @@ module MOM_horizontal_regridding ! character(len=40) :: mdl = "MOM_horizontal_regridding" ! This module's name. !> Fill grid edges -interface fill_boundaries - module procedure fill_boundaries_real - module procedure fill_boundaries_int -end interface +! interface fill_boundaries +! module procedure fill_boundaries_real +! module procedure fill_boundaries_int +! end interface !> Extrapolate and interpolate data interface horiz_interp_and_extrap_tracer @@ -296,7 +292,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, real :: PI_180 integer :: rcode, ncid, varid, ndims, id, jd, kd, jdp - integer :: i,j,k + integer :: i, j, k integer, dimension(4) :: start, count, dims, dim_id real, dimension(:,:), allocatable :: x_in, y_in real, dimension(:), allocatable :: lon_in, lat_in @@ -309,7 +305,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, character(len=8) :: laynum type(horiz_interp_type) :: Interp integer :: is, ie, js, je ! compute domain indices - integer :: isc,iec,jsc,jec ! global compute domain indices + integer :: isc, iec, jsc, jec ! global compute domain indices integer :: isg, ieg, jsg, jeg ! global extent integer :: isd, ied, jsd, jed ! data domain indices integer :: id_clock_read @@ -318,9 +314,9 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, real :: npoints,varAvg real, dimension(SZI_(G),SZJ_(G)) :: lon_out, lat_out, tr_out, mask_out real, dimension(SZI_(G),SZJ_(G)) :: good, fill - real, dimension(SZI_(G),SZJ_(G)) :: tr_outf,tr_prev - real, dimension(SZI_(G),SZJ_(G)) :: good2,fill2 - real, dimension(SZI_(G),SZJ_(G)) :: nlevs + real, dimension(SZI_(G),SZJ_(G)) :: tr_outf, tr_prev + real, dimension(SZI_(G),SZJ_(G)) :: good2, fill2 + real, dimension(SZI_(G),SZJ_(G)) :: nlevs is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed @@ -328,14 +324,14 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, id_clock_read = cpu_clock_id('(Initialize tracer from Z) read', grain=CLOCK_LOOP) - is_ongrid=.false. - if (present(ongrid)) is_ongrid=ongrid + is_ongrid = .false. + if (present(ongrid)) is_ongrid = ongrid if (allocated(tr_z)) deallocate(tr_z) if (allocated(mask_z)) deallocate(mask_z) if (allocated(z_edges_in)) deallocate(z_edges_in) - PI_180=atan(1.0)/45. + PI_180 = atan(1.0)/45. ! Open NetCDF file and if present, extract data and spatial coordinate information ! The convention adopted here requires that the data be written in (i,j,k) ordering. @@ -391,7 +387,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, if (allocated(tr_z)) deallocate(tr_z) if (allocated(mask_z)) deallocate(mask_z) - allocate(lon_in(id),lat_in(jd),z_in(kd),z_edges_in(kd+1)) + allocate(lon_in(id), lat_in(jd), z_in(kd), z_edges_in(kd+1)) allocate(tr_z(isd:ied,jsd:jed,kd), mask_z(isd:ied,jsd:jed,kd)) start = 1 ; count = 1 ; count(1) = id @@ -692,7 +688,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t call cpu_clock_begin(id_clock_read) - fld_sz = get_external_field_size(fms_id) + call get_external_field_info(fms_id, size=fld_sz, axes=axes_data, missing=missing_value) if (allocated(lon_in)) deallocate(lon_in) if (allocated(lat_in)) deallocate(lat_in) if (allocated(z_in)) deallocate(z_in) @@ -700,8 +696,6 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t if (allocated(tr_z)) deallocate(tr_z) if (allocated(mask_z)) deallocate(mask_z) - axes_data = get_external_field_axes(fms_id) - id = fld_sz(1) ; jd = fld_sz(2) ; kd = fld_sz(3) spongeDataOngrid=.false. @@ -722,8 +716,6 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t call cpu_clock_end(id_clock_read) - missing_value = get_external_field_missing(fms_id) - if (.not. spongeDataOngrid) then ! extrapolate the input data to the north pole using the northerm-most latitude max_lat = maxval(lat_in) @@ -773,7 +765,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t if (.not.spongeDataOngrid) then if (is_root_pe()) & - call time_interp_external(fms_id, Time, data_in, verbose=.true., turns=turns) + call time_interp_extern(fms_id, Time, data_in, verbose=.true., turns=turns) ! loop through each data level and interpolate to model grid. ! after interpolating, fill in points which will be needed ! to define the layers @@ -891,7 +883,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t enddo ! kd else - call time_interp_external(fms_id, Time, data_in, verbose=.true., turns=turns) + call time_interp_extern(fms_id, Time, data_in, verbose=.true., turns=turns) do k=1,kd do j=js,je do i=is,ie @@ -927,55 +919,55 @@ end subroutine meshgrid ! None of the subsequent code appears to be used at all. -!> Fill grid edges for integer data -function fill_boundaries_int(m,cyclic_x,tripolar_n) result(mp) - integer, dimension(:,:), intent(in) :: m !< input array (ND) - logical, intent(in) :: cyclic_x !< True if domain is zonally re-entrant - logical, intent(in) :: tripolar_n !< True if domain has an Arctic fold - integer, dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp +! !> Fill grid edges for integer data +! function fill_boundaries_int(m,cyclic_x,tripolar_n) result(mp) +! integer, dimension(:,:), intent(in) :: m !< input array (ND) +! logical, intent(in) :: cyclic_x !< True if domain is zonally re-entrant +! logical, intent(in) :: tripolar_n !< True if domain has an Arctic fold +! integer, dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp - real, dimension(size(m,1),size(m,2)) :: m_real - real, dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp_real +! real, dimension(size(m,1),size(m,2)) :: m_real +! real, dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp_real - m_real = real(m) +! m_real = real(m) - mp_real = fill_boundaries_real(m_real,cyclic_x,tripolar_n) +! mp_real = fill_boundaries_real(m_real,cyclic_x,tripolar_n) - mp = int(mp_real) +! mp = int(mp_real) -end function fill_boundaries_int +! end function fill_boundaries_int !> Fill grid edges for real data -function fill_boundaries_real(m,cyclic_x,tripolar_n) result(mp) - real, dimension(:,:), intent(in) :: m !< input array (ND) - logical, intent(in) :: cyclic_x !< True if domain is zonally re-entrant - logical, intent(in) :: tripolar_n !< True if domain has an Arctic fold - real, dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp +! function fill_boundaries_real(m,cyclic_x,tripolar_n) result(mp) +! real, dimension(:,:), intent(in) :: m !< input array (ND) +! logical, intent(in) :: cyclic_x !< True if domain is zonally re-entrant +! logical, intent(in) :: tripolar_n !< True if domain has an Arctic fold +! real, dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp - integer :: ni,nj,i,j +! integer :: ni,nj,i,j - ni=size(m,1); nj=size(m,2) +! ni=size(m,1); nj=size(m,2) - mp(1:ni,1:nj)=m(:,:) +! mp(1:ni,1:nj)=m(:,:) - if (cyclic_x) then - mp(0,1:nj)=m(ni,1:nj) - mp(ni+1,1:nj)=m(1,1:nj) - else - mp(0,1:nj)=m(1,1:nj) - mp(ni+1,1:nj)=m(ni,1:nj) - endif +! if (cyclic_x) then +! mp(0,1:nj)=m(ni,1:nj) +! mp(ni+1,1:nj)=m(1,1:nj) +! else +! mp(0,1:nj)=m(1,1:nj) +! mp(ni+1,1:nj)=m(ni,1:nj) +! endif - mp(1:ni,0)=m(1:ni,1) - if (tripolar_n) then - do i=1,ni - mp(i,nj+1)=m(ni-i+1,nj) - enddo - else - mp(1:ni,nj+1)=m(1:ni,nj) - endif +! mp(1:ni,0)=m(1:ni,1) +! if (tripolar_n) then +! do i=1,ni +! mp(i,nj+1)=m(ni-i+1,nj) +! enddo +! else +! mp(1:ni,nj+1)=m(1:ni,nj) +! endif -end function fill_boundaries_real +! end function fill_boundaries_real !> Solve del2 (zi) = 0 using successive iterations !! with a 5 point stencil. Only points fill==1 are @@ -983,64 +975,64 @@ end function fill_boundaries_real !! isotropically in index space. The resulting solution !! in each region is an approximation to del2(zi)=0 subject to !! boundary conditions along the valid points curve bounding this region. -subroutine smooth_heights(zi,fill,bad,sor,niter,cyclic_x, tripolar_n) - real, dimension(:,:), intent(inout) :: zi !< input and output array (ND) - integer, dimension(size(zi,1),size(zi,2)), intent(in) :: fill !< same shape as zi, 1=fill - integer, dimension(size(zi,1),size(zi,2)), intent(in) :: bad !< same shape as zi, 1=bad data - real, intent(in) :: sor !< relaxation coefficient (ND) - integer, intent(in) :: niter !< maximum number of iterations - logical, intent(in) :: cyclic_x !< true if domain is zonally reentrant - logical, intent(in) :: tripolar_n !< true if domain has an Arctic fold - - ! Local variables - real, dimension(size(zi,1),size(zi,2)) :: res, m - integer, dimension(size(zi,1),size(zi,2),4) :: B - real, dimension(0:size(zi,1)+1,0:size(zi,2)+1) :: mp - integer, dimension(0:size(zi,1)+1,0:size(zi,2)+1) :: nm - integer :: i,j,k,n - integer :: ni,nj - real :: Isum, bsum - - ni=size(zi,1) ; nj=size(zi,2) - - - mp(:,:) = fill_boundaries(zi,cyclic_x,tripolar_n) - - B(:,:,:) = 0.0 - nm(:,:) = fill_boundaries(bad,cyclic_x,tripolar_n) - - do j=1,nj - do i=1,ni - if (fill(i,j) == 1) then - B(i,j,1)=1-nm(i+1,j);B(i,j,2)=1-nm(i-1,j) - B(i,j,3)=1-nm(i,j+1);B(i,j,4)=1-nm(i,j-1) - endif - enddo - enddo - - do n=1,niter - do j=1,nj - do i=1,ni - if (fill(i,j) == 1) then - bsum = real(B(i,j,1)+B(i,j,2)+B(i,j,3)+B(i,j,4)) - Isum = 1.0/bsum - res(i,j)=Isum*(B(i,j,1)*mp(i+1,j)+B(i,j,2)*mp(i-1,j)+& - B(i,j,3)*mp(i,j+1)+B(i,j,4)*mp(i,j-1)) - mp(i,j) - endif - enddo - enddo - res(:,:)=res(:,:)*sor - - do j=1,nj - do i=1,ni - mp(i,j)=mp(i,j)+res(i,j) - enddo - enddo - - zi(:,:)=mp(1:ni,1:nj) - mp = fill_boundaries(zi,cyclic_x,tripolar_n) - enddo - -end subroutine smooth_heights +! subroutine smooth_heights(zi,fill,bad,sor,niter,cyclic_x, tripolar_n) +! real, dimension(:,:), intent(inout) :: zi !< input and output array (ND) +! integer, dimension(size(zi,1),size(zi,2)), intent(in) :: fill !< same shape as zi, 1=fill +! integer, dimension(size(zi,1),size(zi,2)), intent(in) :: bad !< same shape as zi, 1=bad data +! real, intent(in) :: sor !< relaxation coefficient (ND) +! integer, intent(in) :: niter !< maximum number of iterations +! logical, intent(in) :: cyclic_x !< true if domain is zonally reentrant +! logical, intent(in) :: tripolar_n !< true if domain has an Arctic fold + +! ! Local variables +! real, dimension(size(zi,1),size(zi,2)) :: res, m +! integer, dimension(size(zi,1),size(zi,2),4) :: B +! real, dimension(0:size(zi,1)+1,0:size(zi,2)+1) :: mp +! integer, dimension(0:size(zi,1)+1,0:size(zi,2)+1) :: nm +! integer :: i,j,k,n +! integer :: ni,nj +! real :: Isum, bsum + +! ni=size(zi,1) ; nj=size(zi,2) + + +! mp(:,:) = fill_boundaries(zi,cyclic_x,tripolar_n) + +! B(:,:,:) = 0.0 +! nm(:,:) = fill_boundaries(bad,cyclic_x,tripolar_n) + +! do j=1,nj +! do i=1,ni +! if (fill(i,j) == 1) then +! B(i,j,1)=1-nm(i+1,j);B(i,j,2)=1-nm(i-1,j) +! B(i,j,3)=1-nm(i,j+1);B(i,j,4)=1-nm(i,j-1) +! endif +! enddo +! enddo + +! do n=1,niter +! do j=1,nj +! do i=1,ni +! if (fill(i,j) == 1) then +! bsum = real(B(i,j,1)+B(i,j,2)+B(i,j,3)+B(i,j,4)) +! Isum = 1.0/bsum +! res(i,j)=Isum*(B(i,j,1)*mp(i+1,j)+B(i,j,2)*mp(i-1,j)+& +! B(i,j,3)*mp(i,j+1)+B(i,j,4)*mp(i,j-1)) - mp(i,j) +! endif +! enddo +! enddo +! res(:,:)=res(:,:)*sor + +! do j=1,nj +! do i=1,ni +! mp(i,j)=mp(i,j)+res(i,j) +! enddo +! enddo + +! zi(:,:)=mp(1:ni,1:nj) +! mp = fill_boundaries(zi,cyclic_x,tripolar_n) +! enddo + +! end subroutine smooth_heights end module MOM_horizontal_regridding diff --git a/src/framework/MOM_interpolate.F90 b/src/framework/MOM_interpolate.F90 new file mode 100644 index 0000000000..c63c847e55 --- /dev/null +++ b/src/framework/MOM_interpolate.F90 @@ -0,0 +1,143 @@ +!> This module wraps the FMS temporal and spatial interpolation routines +module MOM_interpolate + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_array_transform, only : allocate_rotated_array, rotate_array +use MOM_error_handler, only : MOM_error, FATAL +use MOM_io_wrapper, only : axistype +use horiz_interp_mod, only : horiz_interp_new, horiz_interp, horiz_interp_init, horiz_interp_type +use time_interp_external_mod, only : time_interp_external_fms=>time_interp_external +use time_interp_external_mod, only : init_external_field, time_interp_external_init +use time_interp_external_mod, only : get_external_field_size +use time_interp_external_mod, only : get_external_field_axes, get_external_field_missing +use time_manager_mod, only : time_type + +implicit none ; private + +public :: time_interp_extern, init_external_field, time_interp_external_init +public :: get_external_field_info +public :: horiz_interp_type, horiz_interp_init, horiz_interp, horiz_interp_new + +!> Read a field based on model time, and rotate to the model domain. +! This inerface does not share the name time_interp_external with the module it primarily +! wraps because of errors (perhaps a bug) that arise with the PGI 19.10.0 compiler. +interface time_interp_extern + module procedure time_interp_external_0d + module procedure time_interp_external_2d + module procedure time_interp_external_3d +end interface time_interp_extern + +contains + +!> Get information about the external fields. +subroutine get_external_field_info(field_id, size, axes, missing) + integer, intent(in) :: field_id !< The integer index of the external + !! field returned from a previous + !! call to init_external_field() + integer, dimension(4), optional, intent(inout) :: size !< Dimension sizes for the input data + type(axistype), dimension(4), optional, intent(inout) :: axes !< Axis types for the input data + real, optional, intent(inout) :: missing !< Missing value for the input data + + if (present(size)) then + size(1:4) = get_external_field_size(field_id) + endif + + if (present(axes)) then + axes(1:4) = get_external_field_axes(field_id) + endif + + if (present(missing)) then + missing = get_external_field_missing(field_id) + endif + +end subroutine get_external_field_info + + +!> Read a scalar field based on model time. +subroutine time_interp_external_0d(field_id, time, data_in, verbose) + integer, intent(in) :: field_id !< The integer index of the external field returned + !! from a previous call to init_external_field() + type(time_type), intent(in) :: time !< The target time for the data + real, intent(inout) :: data_in !< The interpolated value + logical, optional, intent(in) :: verbose !< If true, write verbose output for debugging + + call time_interp_external_fms(field_id, time, data_in, verbose=verbose) +end subroutine time_interp_external_0d + +!> Read a 2d field from an external based on model time, potentially including horizontal +!! interpolation and rotation of the data +subroutine time_interp_external_2d(field_id, time, data_in, interp, verbose, horz_interp, mask_out, turns) + integer, intent(in) :: field_id !< The integer index of the external field returned + !! from a previous call to init_external_field() + type(time_type), intent(in) :: time !< The target time for the data + real, dimension(:,:), intent(inout) :: data_in !< The array in which to store the interpolated values + integer, optional, intent(in) :: interp !< A flag indicating the temporal interpolation method + logical, optional, intent(in) :: verbose !< If true, write verbose output for debugging + type(horiz_interp_type), & + optional, intent(in) :: horz_interp !< A structure to control horizontal interpolation + logical, dimension(:,:), & + optional, intent(out) :: mask_out !< An array that is true where there is valid data + integer, optional, intent(in) :: turns !< Number of quarter turns to rotate the data + + real, allocatable :: data_pre_rot(:,:) ! The input data before rotation + integer :: qturns ! The number of quarter turns to rotate the data + + ! TODO: Mask rotation requires logical array rotation support + if (present(mask_out)) & + call MOM_error(FATAL, "Rotation of masked output not yet support") + + qturns = 0 + if (present(turns)) qturns = modulo(turns, 4) + + if (qturns == 0) then + call time_interp_external_fms(field_id, time, data_in, interp=interp, & + verbose=verbose, horz_interp=horz_interp) + else + call allocate_rotated_array(data_in, [1,1], -qturns, data_pre_rot) + call time_interp_external_fms(field_id, time, data_pre_rot, interp=interp, & + verbose=verbose, horz_interp=horz_interp) + call rotate_array(data_pre_rot, turns, data_in) + deallocate(data_pre_rot) + endif +end subroutine time_interp_external_2d + + +!> Read a 3d field based on model time, and rotate to the model grid +subroutine time_interp_external_3d(field_id, time, data_in, interp, & + verbose, horz_interp, mask_out, turns) + integer, intent(in) :: field_id !< The integer index of the external field returned + !! from a previous call to init_external_field() + type(time_type), intent(in) :: time !< The target time for the data + real, dimension(:,:,:), intent(inout) :: data_in !< The array in which to store the interpolated values + integer, optional, intent(in) :: interp !< A flag indicating the temporal interpolation method + logical, optional, intent(in) :: verbose !< If true, write verbose output for debugging + type(horiz_interp_type), & + optional, intent(in) :: horz_interp !< A structure to control horizontal interpolation + logical, dimension(:,:,:), & + optional, intent(out) :: mask_out !< An array that is true where there is valid data + integer, optional, intent(in) :: turns !< Number of quarter turns to rotate the data + + real, allocatable :: data_pre_rot(:,:,:) ! The input data before rotation + integer :: qturns ! The number of quarter turns to rotate the data + + ! TODO: Mask rotation requires logical array rotation support + if (present(mask_out)) & + call MOM_error(FATAL, "Rotation of masked output not yet support") + + qturns = 0 + if (present(turns)) qturns = modulo(turns, 4) + + if (qturns == 0) then + call time_interp_external_fms(field_id, time, data_in, interp=interp, & + verbose=verbose, horz_interp=horz_interp) + else + call allocate_rotated_array(data_in, [1,1,1], -qturns, data_pre_rot) + call time_interp_external_fms(field_id, time, data_pre_rot, interp=interp, & + verbose=verbose, horz_interp=horz_interp) + call rotate_array(data_pre_rot, turns, data_in) + deallocate(data_pre_rot) + endif +end subroutine time_interp_external_3d + +end module MOM_interpolate diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index a5dd92b640..634ff80325 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -57,8 +57,8 @@ module MOM_ice_shelf use MOM_coms, only : reproducing_sum use MOM_spatial_means, only : global_area_integral use MOM_checksums, only : hchksum, qchksum, chksum, uchksum, vchksum, uvchksum -use time_interp_external_mod, only : init_external_field, time_interp_external -use time_interp_external_mod, only : time_interp_external_init +use MOM_interpolate, only : init_external_field, time_interp_extern, time_interp_external_init + implicit none ; private #include @@ -1084,9 +1084,9 @@ subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes) do j=js,je ; do i=is,ie last_hmask(i,j) = ISS%hmask(i,j) ; last_area_shelf_h(i,j) = ISS%area_shelf_h(i,j) enddo ; enddo - call time_interp_external(CS%id_read_mass, Time0, last_mass_shelf) + call time_interp_extern(CS%id_read_mass, Time0, last_mass_shelf) do j=js,je ; do i=is,ie - ! This should only be done if time_interp_external did an update. + ! This should only be done if time_interp_extern did an update. last_mass_shelf(i,j) = US%kg_m3_to_R*US%m_to_Z * last_mass_shelf(i,j) ! Rescale after time_interp last_h_shelf(i,j) = last_mass_shelf(i,j) / CS%density_ice enddo ; enddo @@ -1512,7 +1512,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, if (CS%rotate_index) then allocate(tmp2d(CS%Grid_in%isd:CS%Grid_in%ied,CS%Grid_in%jsd:CS%Grid_in%jed));tmp2d(:,:)=0.0 call MOM_read_data(TideAmp_file, 'tideamp', tmp2d, CS%Grid_in%domain, timelevel=1, scale=US%m_s_to_L_T) - call rotate_array(tmp2d,CS%turns, CS%utide) + call rotate_array(tmp2d, CS%turns, CS%utide) deallocate(tmp2d) else call MOM_read_data(TideAmp_file, 'tideamp', CS%utide, CS%Grid%domain, timelevel=1, scale=US%m_s_to_L_T) @@ -1984,8 +1984,8 @@ subroutine update_shelf_mass(G, US, CS, ISS, Time) allocate(tmp2d(is:ie,js:je)) ; tmp2d(:,:) = 0.0 endif - call time_interp_external(CS%id_read_mass, Time, tmp2d) - call rotate_array(tmp2d,CS%turns, ISS%mass_shelf) + call time_interp_extern(CS%id_read_mass, Time, tmp2d) + call rotate_array(tmp2d, CS%turns, ISS%mass_shelf) deallocate(tmp2d) ! This should only be done if time_interp_external did an update. diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 470098a08a..8ba9dd959a 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -15,15 +15,15 @@ module MOM_diabatic_aux use MOM_file_parser, only : get_param, log_param, log_version, param_file_type use MOM_forcing_type, only : forcing, extractFluxes1d, forcing_SinglePointPrint use MOM_grid, only : ocean_grid_type +use MOM_interpolate, only : init_external_field, time_interp_extern +use MOM_interpolate, only : time_interp_external_init use MOM_io, only : slasher use MOM_opacity, only : set_opacity, opacity_CS, extract_optics_slice, extract_optics_fields use MOM_opacity, only : optics_type, optics_nbands, absorbRemainingSW, sumSWoverBands use MOM_tracer_flow_control, only : get_chl_from_model, tracer_flow_control_CS use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : thermo_var_ptrs ! , vertvisc_type, accel_diag_ptrs +use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type -use time_interp_external_mod, only : init_external_field, time_interp_external -use time_interp_external_mod, only : time_interp_external_init implicit none ; private @@ -621,7 +621,7 @@ subroutine set_pen_shortwave(optics, fluxes, G, GV, US, CS, opacity_CSp, tracer_ if (CS%chl_from_file) then ! Only the 2-d surface chlorophyll can be read in from a file. The ! same value is assumed for all layers. - call time_interp_external(CS%sbc_chl, CS%Time, chl_2d) + call time_interp_extern(CS%sbc_chl, CS%Time, chl_2d) do j=js,je ; do i=is,ie if ((G%mask2dT(i,j) > 0.5) .and. (chl_2d(i,j) < 0.0)) then write(mesg,'(" Time_interp negative chl of ",(1pe12.4)," at i,j = ",& From 43ae9ae7d6340db6a2500acdaf0d7410de2c3462 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Wed, 13 Jan 2021 18:29:16 -0500 Subject: [PATCH 05/11] Adds --with-driver option to configure - Following suggestion from @marshallward I've implemented an option to select the driver code at the configure stage (changes to configure.ac) - This allows the target build/coupled/ocean_model_MOM.o to be built using the .testing ac tools. --- .testing/Makefile | 23 ++++++----------------- ac/configure.ac | 10 ++++++++-- 2 files changed, 14 insertions(+), 19 deletions(-) diff --git a/.testing/Makefile b/.testing/Makefile index a780f19323..2806d54130 100644 --- a/.testing/Makefile +++ b/.testing/Makefile @@ -226,7 +226,9 @@ build/asymmetric/Makefile: MOM_ENV=$(PATH_FMS) $(ASYMMETRIC_FCFLAGS) $(MOM_LDFLA build/repro/Makefile: MOM_ENV=$(PATH_FMS) $(REPRO_FCFLAGS) $(MOM_LDFLAGS) build/openmp/Makefile: MOM_ENV=$(PATH_FMS) $(OPENMP_FCFLAGS) $(MOM_LDFLAGS) build/target/Makefile: MOM_ENV=$(PATH_FMS) $(TARGET_FCFLAGS) $(MOM_LDFLAGS) - +build/coupled/Makefile: MOM_ENV=$(PATH_FMS) $(SYMMETRIC_FCFLAGS) $(SYMMETRIC_LDFLAGS) +build/nuopc/Makefile: MOM_ENV=$(PATH_FMS) $(SYMMETRIC_FCFLAGS) $(SYMMETRIC_LDFLAGS) +build/mct/Makefile: MOM_ENV=$(PATH_FMS) $(SYMMETRIC_FCFLAGS) $(SYMMETRIC_LDFLAGS) # Configure script flags build/symmetric/Makefile: MOM_ACFLAGS= @@ -234,7 +236,9 @@ build/asymmetric/Makefile: MOM_ACFLAGS=--enable-asymmetric build/repro/Makefile: MOM_ACFLAGS= build/openmp/Makefile: MOM_ACFLAGS=--enable-openmp build/target/Makefile: MOM_ACFLAGS= - +build/coupled/Makefile: MOM_ACFLAGS=--with-driver=coupled_driver +build/nuopc/Makefile: MOM_ACFLAGS=--with-driver=nuopc_driver +build/mct/Makefile: MOM_ACFLAGS=--with-driver=mct_driver # Fetch regression target source code build/target/Makefile: | $(TARGET_CODEBASE) @@ -336,29 +340,14 @@ $(DEPS)/Makefile: ../ac/deps/Makefile # - use autoconf rather than mkmf templates MK_TEMPLATE ?= ../../$(DEPS)/mkmf/templates/ncrc-gnu.mk # NUOPC driver -build/nuopc/Makefile: ../config_src/nuopc_driver $(MOM_SOURCE) - mkdir -p $(@D) - cd $(@D) \ - && $(MOM_ENV) ../../$(DEPS)/mkmf/bin/list_paths -l ../../$(DEPS)/fms/src ../../../config_src/nuopc_driver ../../../config_src/dynamic_symmetric ../../../src/ ../../../config_src/external \ - && $(MOM_ENV) ../../$(DEPS)/mkmf/bin/mkmf -t $(MK_TEMPLATE) -p MOM6 path_names build/nuopc/mom_ocean_model_nuopc.o: build/nuopc/Makefile cd $(@D) && make $(@F) check_mom6_api_nuopc: build/nuopc/mom_ocean_model_nuopc.o # GFDL coupled driver -build/coupled/Makefile: ../config_src/coupled_driver $(MOM_SOURCE) - mkdir -p $(@D) - cd $(@D) \ - && $(MOM_ENV) ../../$(DEPS)/mkmf/bin/list_paths -l ../../$(DEPS)/fms/src ../../../config_src/coupled_driver ../../../config_src/dynamic_symmetric ../../../src/ ../../../config_src/external \ - && $(MOM_ENV) ../../$(DEPS)/mkmf/bin/mkmf -t $(MK_TEMPLATE) -p MOM6 path_names build/coupled/ocean_model_MOM.o: build/coupled/Makefile cd $(@D) && make $(@F) check_mom6_api_coupled: build/coupled/ocean_model_MOM.o # MCT driver -build/mct/Makefile: ../config_src/mct_driver $(MOM_SOURCE) - mkdir -p $(@D) - cd $(@D) \ - && $(MOM_ENV) ../../$(DEPS)/mkmf/bin/list_paths -l ../../$(DEPS)/fms/src ../../../config_src/mct_driver ../../../config_src/dynamic_symmetric ../../../src/ ../../../config_src/external \ - && $(MOM_ENV) ../../$(DEPS)/mkmf/bin/mkmf -t $(MK_TEMPLATE) -p MOM6 path_names build/mct/mom_ocean_model_mct.o: build/mct/Makefile cd $(@D) && make $(@F) check_mom6_api_mct: build/mct/mom_ocean_model_mct.o diff --git a/ac/configure.ac b/ac/configure.ac index ad8ed83603..487230beb8 100644 --- a/ac/configure.ac +++ b/ac/configure.ac @@ -48,6 +48,12 @@ AC_ARG_ENABLE([asymmetric], AS_IF([test "$enable_asymmetric" = yes], [MEM_LAYOUT=${srcdir}/config_src/dynamic]) +# Default to solo_driver +DRIVER_DIR=${srcdir}/config_src/solo_driver +AC_ARG_WITH([driver], + AS_HELP_STRING([--with-driver=coupled_driver|solo_driver], [Select directory for driver source code])) +AS_IF([test "x$with_driver" != "x"], + [DRIVER_DIR=${srcdir}/config_src/${with_driver}]) # TODO: Rather than point to a pre-configured header file, autoconf could be # used to configure a header based on a template. @@ -210,10 +216,10 @@ AS_IF([test -z "$MKMF"], [ AC_CONFIG_COMMANDS([path_names], [list_paths -l \ ${srcdir}/src \ - ${srcdir}/config_src/solo_driver \ ${srcdir}/config_src/ext* \ + ${DRIVER_DIR} \ ${MEM_LAYOUT} -], [MEM_LAYOUT=$MEM_LAYOUT]) +], [MEM_LAYOUT=$MEM_LAYOUT DRIVER_DIR=$DRIVER_DIR]) AC_CONFIG_COMMANDS([Makefile.mkmf], From 76cb4718d01dab9ab0c714cba310f13675c5f5dc Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 13 Jan 2021 16:03:32 -0500 Subject: [PATCH 06/11] Cleaned up MOM_horizontal_regridding Partially cleaned up MOM_horizontal_regridding.F90, including the elimination of several unused routines, the addition of comments describing variables, and some white space corrections. All answers are bitwise identical. --- src/framework/MOM_horizontal_regridding.F90 | 335 +++++++------------- 1 file changed, 110 insertions(+), 225 deletions(-) diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index 537c7e4d92..f1ab073938 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -11,10 +11,10 @@ module MOM_horizontal_regridding use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint use MOM_file_parser, only : get_param, log_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type -use MOM_io_wrapper, only : axistype, get_axis_data -use MOM_time_manager, only : time_type use MOM_interpolate, only : time_interp_extern, get_external_field_info, horiz_interp_init use MOM_interpolate, only : horiz_interp_new, horiz_interp, horiz_interp_type +use MOM_io_wrapper, only : axistype, get_axis_data +use MOM_time_manager, only : time_type use netcdf @@ -24,14 +24,6 @@ module MOM_horizontal_regridding public :: horiz_interp_and_extrap_tracer, myStats -! character(len=40) :: mdl = "MOM_horizontal_regridding" ! This module's name. - -!> Fill grid edges -! interface fill_boundaries -! module procedure fill_boundaries_real -! module procedure fill_boundaries_int -! end interface - !> Extrapolate and interpolate data interface horiz_interp_and_extrap_tracer module procedure horiz_interp_and_extrap_tracer_record @@ -80,10 +72,9 @@ subroutine myStats(array, missing, is, ie, js, je, k, mesg) end subroutine myStats -!> Use ICE-9 algorithm to populate points (fill=1) with -!! valid data (good=1). If no information is available, -!! Then use a previous guess (prev). Optionally (smooth) -!! blend the filled points to achieve a more desirable result. +!> Use ICE-9 algorithm to populate points (fill=1) with valid data (good=1). If no information +!! is available, use a previous guess (prev). Optionally (smooth) blend the filled points to +!! achieve a more desirable result. subroutine fill_miss_2d(aout, good, fill, prev, G, smooth, num_pass, relc, crit, debug, answers_2018) use MOM_coms, only : sum_across_PEs @@ -108,15 +99,20 @@ subroutine fill_miss_2d(aout, good, fill, prev, G, smooth, num_pass, relc, crit, !! answers as the code did in late 2018. Otherwise !! add parentheses for rotational symmetry. - real, dimension(SZI_(G),SZJ_(G)) :: b,r - real, dimension(SZI_(G),SZJ_(G)) :: fill_pts, good_, good_new - + real, dimension(SZI_(G),SZJ_(G)) :: a_filled ! The aout with missing values filled in + real, dimension(SZI_(G),SZJ_(G)) :: a_chg ! The change in aout due to an iteration of smoothing + real, dimension(SZI_(G),SZJ_(G)) :: fill_pts ! 1 for points that still need to be filled + real, dimension(SZI_(G),SZJ_(G)) :: good_ ! The values that are valid for the current iteration + real, dimension(SZI_(G),SZJ_(G)) :: good_new ! The values of good_ to use for the next iteration + + real :: east, west, north, south ! Valid neighboring values or 0 for invalid values + real :: ge, gw, gn, gs ! Flags indicating which neighbors have valid values + real :: ngood ! The number of valid values in neighboring points + logical :: do_smooth ! Indicates whether to do smoothing of the array + real :: nfill ! The remaining number of points to fill + real :: nfill_prev ! The previous value of nfill character(len=256) :: mesg ! The text of an error message - integer :: i,j,k - real :: east,west,north,south,sor - real :: ge,gw,gn,gs,ngood - logical :: do_smooth,siena_bug - real :: nfill, nfill_prev + integer :: i, j, k integer, parameter :: num_pass_default = 10000 real, parameter :: relc_default = 0.25, crit_default = 1.e-3 @@ -151,15 +147,15 @@ subroutine fill_miss_2d(aout, good, fill, prev, G, smooth, num_pass, relc, crit, nfill_prev = nfill good_(:,:) = good(:,:) - r(:,:) = 0.0 + a_chg(:,:) = 0.0 do while (nfill > 0.0) call pass_var(good_,G%Domain) call pass_var(aout,G%Domain) - b(:,:)=aout(:,:) - good_new(:,:)=good_(:,:) + a_filled(:,:) = aout(:,:) + good_new(:,:) = good_(:,:) do j=js,je ; do i=is,ie @@ -180,16 +176,16 @@ subroutine fill_miss_2d(aout, good, fill, prev, G, smooth, num_pass, relc, crit, endif if (ngood > 0.) then if (ans_2018) then - b(i,j)=(east+west+north+south)/ngood + a_filled(i,j) = (east+west+north+south)/ngood else - b(i,j) = ((east+west) + (north+south))/ngood + a_filled(i,j) = ((east+west) + (north+south))/ngood endif fill_pts(i,j) = 0.0 good_new(i,j) = 1.0 endif enddo ; enddo - aout(is:ie,js:je) = b(is:ie,js:je) + aout(is:ie,js:je) = a_filled(is:ie,js:je) good_(is:ie,js:je) = good_new(is:ie,js:je) nfill_prev = nfill nfill = sum(fill_pts(is:ie,js:je)) @@ -209,11 +205,13 @@ subroutine fill_miss_2d(aout, good, fill, prev, G, smooth, num_pass, relc, crit, call MOM_error(WARNING, mesg, .true.) endif + ! Determine the number of remaining points to fill globally. nfill = sum(fill_pts(is:ie,js:je)) call sum_across_PEs(nfill) - enddo + enddo ! while block for remaining points to fill. + ! Do Laplacian smoothing for the points that have been filled in. if (do_smooth) then ; do k=1,npass call pass_var(aout,G%Domain) do j=js,je ; do i=is,ie @@ -221,22 +219,22 @@ subroutine fill_miss_2d(aout, good, fill, prev, G, smooth, num_pass, relc, crit, east = max(good(i+1,j),fill(i+1,j)) ; west = max(good(i-1,j),fill(i-1,j)) north = max(good(i,j+1),fill(i,j+1)) ; south = max(good(i,j-1),fill(i,j-1)) if (ans_2018) then - r(i,j) = relax_coeff*(south*aout(i,j-1)+north*aout(i,j+1) + & - west*aout(i-1,j)+east*aout(i+1,j) - & - (south+north+west+east)*aout(i,j)) + a_chg(i,j) = relax_coeff*(south*aout(i,j-1)+north*aout(i,j+1) + & + west*aout(i-1,j)+east*aout(i+1,j) - & + (south+north+west+east)*aout(i,j)) else - r(i,j) = relax_coeff*( ((south*aout(i,j-1) + north*aout(i,j+1)) + & + a_chg(i,j) = relax_coeff*( ((south*aout(i,j-1) + north*aout(i,j+1)) + & (west*aout(i-1,j)+east*aout(i+1,j))) - & ((south+north)+(west+east))*aout(i,j) ) endif else - r(i,j) = 0. + a_chg(i,j) = 0. endif enddo ; enddo ares = 0.0 do j=js,je ; do i=is,ie - aout(i,j) = r(i,j) + aout(i,j) - ares = max(ares, abs(r(i,j))) + aout(i,j) = a_chg(i,j) + aout(i,j) + ares = max(ares, abs(a_chg(i,j))) enddo ; enddo call max_across_PEs(ares) if (ares <= acrit) exit @@ -285,9 +283,9 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, !! extrapolation is performed by this routine ! Local variables - real, dimension(:,:), allocatable :: tr_in, tr_inp ! A 2-d array for holding input data on - ! native horizontal grid and extended grid - ! with poles. + real, dimension(:,:), allocatable :: tr_in !< A 2-d array for holding input data on its + !! native horizontal grid. + real, dimension(:,:), allocatable :: tr_inp !< Native horizontal grid data extended to the poles. real, dimension(:,:), allocatable :: mask_in ! A 2-d mask for extended input grid. real :: PI_180 @@ -295,8 +293,8 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, integer :: i, j, k integer, dimension(4) :: start, count, dims, dim_id real, dimension(:,:), allocatable :: x_in, y_in - real, dimension(:), allocatable :: lon_in, lat_in - real, dimension(:), allocatable :: lat_inp, last_row + real, dimension(:), allocatable :: lon_in, lat_in ! The longitude and latitude in the input file + real, dimension(:), allocatable :: lat_inp ! The input file latitudes expanded to the pole real :: max_lat, min_lat, pole, max_depth, npole real :: roundoff ! The magnitude of roundoff, usually ~2e-16. real :: add_offset, scale_factor @@ -311,12 +309,15 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, integer :: id_clock_read character(len=12) :: dim_name(4) logical :: debug=.false. - real :: npoints,varAvg - real, dimension(SZI_(G),SZJ_(G)) :: lon_out, lat_out, tr_out, mask_out - real, dimension(SZI_(G),SZJ_(G)) :: good, fill - real, dimension(SZI_(G),SZJ_(G)) :: tr_outf, tr_prev - real, dimension(SZI_(G),SZJ_(G)) :: good2, fill2 - real, dimension(SZI_(G),SZJ_(G)) :: nlevs + real :: npoints, varAvg + real, dimension(SZI_(G),SZJ_(G)) :: lon_out, lat_out ! The longitude and latitude of points on the model grid + real, dimension(SZI_(G),SZJ_(G)) :: tr_out, mask_out ! The tracer and mask on the model grid + real, dimension(SZI_(G),SZJ_(G)) :: good ! Where the data is valid, this is 1. + real, dimension(SZI_(G),SZJ_(G)) :: fill ! 1 where the data needs to be filled in + real, dimension(SZI_(G),SZJ_(G)) :: tr_outf ! The tracer concentrations after Ice-9 + real, dimension(SZI_(G),SZJ_(G)) :: tr_prev ! The tracer concentrations in the layer above + real, dimension(SZI_(G),SZJ_(G)) :: good2 ! 1 where the data is valid after Ice-9 + real, dimension(SZI_(G),SZJ_(G)) :: fill2 ! 1 for points that still need to be filled after Ice-9 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed @@ -407,7 +408,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, if (present(m_to_Z)) then ; do k=1,kd ; z_in(k) = m_to_Z * z_in(k) ; enddo ; endif - ! extrapolate the input data to the north pole using the northerm-most latitude + ! extrapolate the input data to the north pole using the northern-most latitude add_np = .false. jdp = jd if (.not. is_ongrid) then @@ -445,7 +446,6 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, allocate(tr_in(id,jd)) ; tr_in(:,:) = 0.0 allocate(tr_inp(id,jdp)) ; tr_inp(:,:) = 0.0 allocate(mask_in(id,jdp)) ; mask_in(:,:) = 0.0 - allocate(last_row(id)) ; last_row(:) = 0.0 endif max_depth = maxval(G%bathyT) @@ -488,11 +488,11 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, trim(varnam)//" in file "// trim(filename)) if (add_np) then - last_row(:)=tr_in(:,jd); pole=0.0;npole=0.0 + pole = 0.0 ; npole = 0.0 do i=1,id if (abs(tr_in(i,jd)-missing_value) > abs(roundoff*missing_value)) then - pole = pole+last_row(i) - npole = npole+1.0 + pole = pole + tr_in(i,jd) + npole = npole + 1.0 endif enddo if (npole > 0) then @@ -525,12 +525,12 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, ! call fms routine horiz_interp to interpolate input level data to model horizontal grid if (.not. is_ongrid) then if (k == 1) then - call horiz_interp_new(Interp,x_in,y_in,lon_out(is:ie,js:je),lat_out(is:ie,js:je), & - interp_method='bilinear',src_modulo=.true.) + call horiz_interp_new(Interp, x_in, y_in, lon_out(is:ie,js:je), lat_out(is:ie,js:je), & + interp_method='bilinear', src_modulo=.true.) endif if (debug) then - call myStats(tr_inp,missing_value, is,ie,js,je,k,'Tracer from file') + call myStats(tr_inp,missing_value, is, ie, js, je, k,'Tracer from file') endif endif @@ -538,7 +538,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, if (is_ongrid) then tr_out(is:ie,js:je)=tr_in(is:ie,js:je) else - call horiz_interp(Interp,tr_inp,tr_out(is:ie,js:je), missing_value=missing_value, new_missing_handle=.true.) + call horiz_interp(Interp, tr_inp, tr_out(is:ie,js:je), missing_value=missing_value, new_missing_handle=.true.) endif mask_out=1.0 @@ -591,7 +591,9 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, fill2(:,:) = fill(:,:) call fill_miss_2d(tr_outf, good2, fill2, tr_prev, G, smooth=.true., answers_2018=answers_2018) - call myStats(tr_outf, missing_value, is, ie, js, je, k, 'field from fill_miss_2d()') + if (debug) then + call myStats(tr_outf, missing_value, is, ie, js, je, k, 'field from fill_miss_2d()') + endif tr_z(:,:,k) = tr_outf(:,:) * G%mask2dT(:,:) mask_z(:,:,k) = good2(:,:) + fill2(:,:) @@ -634,9 +636,9 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t !! add parentheses for rotational symmetry. ! Local variables - real, dimension(:,:), allocatable :: tr_in,tr_inp !< A 2-d array for holding input data on - !! native horizontal grid and extended grid - !! with poles. + real, dimension(:,:), allocatable :: tr_in !< A 2-d array for holding input data on its + !! native horizontal grid. + real, dimension(:,:), allocatable :: tr_inp !< Native horizontal grid data extended to the poles. real, dimension(:,:,:), allocatable :: data_in !< A buffer for storing the full 3-d time-interpolated array !! on the original grid real, dimension(:,:), allocatable :: mask_in !< A 2-d mask for extended input grid. @@ -646,8 +648,8 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t integer :: i,j,k integer, dimension(4) :: start, count, dims, dim_id real, dimension(:,:), allocatable :: x_in, y_in - real, dimension(:), allocatable :: lon_in, lat_in - real, dimension(:), allocatable :: lat_inp, last_row + real, dimension(:), allocatable :: lon_in, lat_in ! The longitude and latitude in the input file + real, dimension(:), allocatable :: lat_inp ! The input file latitudes expanded to the pole real :: max_lat, min_lat, pole, max_depth, npole real :: roundoff ! The magnitude of roundoff, usually ~2e-16. logical :: add_np @@ -663,12 +665,15 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t character(len=12) :: dim_name(4) logical :: debug=.false. logical :: spongeDataOngrid - real :: npoints,varAvg - real, dimension(SZI_(G),SZJ_(G)) :: lon_out, lat_out, tr_out, mask_out - real, dimension(SZI_(G),SZJ_(G)) :: good, fill - real, dimension(SZI_(G),SZJ_(G)) :: tr_outf,tr_prev - real, dimension(SZI_(G),SZJ_(G)) :: good2,fill2 - real, dimension(SZI_(G),SZJ_(G)) :: nlevs + real :: npoints, varAvg + real, dimension(SZI_(G),SZJ_(G)) :: lon_out, lat_out ! The longitude and latitude of points on the model grid + real, dimension(SZI_(G),SZJ_(G)) :: tr_out, mask_out ! The tracer and mask on the model grid + real, dimension(SZI_(G),SZJ_(G)) :: good ! Where the data is valid, this is 1. + real, dimension(SZI_(G),SZJ_(G)) :: fill ! 1 where the data needs to be filled in + real, dimension(SZI_(G),SZJ_(G)) :: tr_outf ! The tracer concentrations after Ice-9 + real, dimension(SZI_(G),SZJ_(G)) :: tr_prev ! The tracer concentrations in the layer above + real, dimension(SZI_(G),SZJ_(G)) :: good2 ! 1 where the data is valid after Ice-9 + real, dimension(SZI_(G),SZJ_(G)) :: fill2 ! 1 for points that still need to be filled after Ice-9 integer :: turns turns = G%HI%turns @@ -679,7 +684,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t id_clock_read = cpu_clock_id('(Initialize tracer from Z) read', grain=CLOCK_LOOP) - PI_180=atan(1.0)/45. + PI_180 = atan(1.0)/45. ! Open NetCDF file and if present, extract data and spatial coordinate information ! The convention adopted here requires that the data be written in (i,j,k) ordering. @@ -696,15 +701,15 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t id = fld_sz(1) ; jd = fld_sz(2) ; kd = fld_sz(3) - spongeDataOngrid=.false. - if (PRESENT(spongeOngrid)) spongeDataOngrid=spongeOngrid + spongeDataOngrid = .false. + if (PRESENT(spongeOngrid)) spongeDataOngrid = spongeOngrid if (.not. spongeDataOngrid) then - allocate(lon_in(id),lat_in(jd)) + allocate(lon_in(id), lat_in(jd)) call get_axis_data(axes_data(1), lon_in) call get_axis_data(axes_data(2), lat_in) endif - allocate(z_in(kd),z_edges_in(kd+1)) + allocate(z_in(kd), z_edges_in(kd+1)) allocate(tr_z(isd:ied,jsd:jed,kd), mask_z(isd:ied,jsd:jed,kd)) @@ -715,9 +720,9 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t call cpu_clock_end(id_clock_read) if (.not. spongeDataOngrid) then - ! extrapolate the input data to the north pole using the northerm-most latitude + ! Extrapolate the input data to the north pole using the northerm-most latitude. max_lat = maxval(lat_in) - add_np=.false. + add_np = .false. if (max_lat < 90.0) then add_np = .true. jdp = jd+1 @@ -728,7 +733,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t allocate(lat_in(1:jdp)) lat_in(:) = lat_inp(:) else - jdp=jd + jdp = jd endif call horiz_interp_init() lon_in = lon_in*PI_180 @@ -741,7 +746,6 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t allocate(tr_in(id,jd)) ; tr_in(:,:)=0.0 allocate(tr_inp(id,jdp)) ; tr_inp(:,:)=0.0 allocate(mask_in(id,jdp)) ; mask_in(:,:)=0.0 - allocate(last_row(id)) ; last_row(:)=0.0 else allocate(data_in(isd:ied,jsd:jed,kd)) endif @@ -764,40 +768,39 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t if (.not.spongeDataOngrid) then if (is_root_pe()) & call time_interp_extern(fms_id, Time, data_in, verbose=.true., turns=turns) - ! loop through each data level and interpolate to model grid. - ! after interpolating, fill in points which will be needed - ! to define the layers + ! Loop through each data level and interpolate to model grid. + ! After interpolating, fill in points which will be needed to define the layers. do k=1,kd write(laynum,'(I8)') k ; laynum = adjustl(laynum) if (is_root_pe()) then tr_in(1:id,1:jd) = data_in(1:id,1:jd,k) if (add_np) then - last_row(:)=tr_in(:,jd); pole=0.0;npole=0.0 - do i=1,id - if (abs(tr_in(i,jd)-missing_value) > abs(roundoff*missing_value)) then - pole = pole+last_row(i) - npole = npole+1.0 - endif - enddo - if (npole > 0) then - pole=pole/npole - else - pole=missing_value - endif - tr_inp(:,1:jd) = tr_in(:,:) - tr_inp(:,jdp) = pole + pole = 0.0 ; npole = 0.0 + do i=1,id + if (abs(tr_in(i,jd)-missing_value) > abs(roundoff*missing_value)) then + pole = pole + tr_in(i,jd) + npole = npole+1.0 + endif + enddo + if (npole > 0) then + pole = pole / npole + else + pole = missing_value + endif + tr_inp(:,1:jd) = tr_in(:,:) + tr_inp(:,jdp) = pole else - tr_inp(:,:) = tr_in(:,:) + tr_inp(:,:) = tr_in(:,:) endif endif call broadcast(tr_inp, id*jdp, blocking=.true.) - mask_in=0.0 + mask_in(:,:) = 0.0 do j=1,jdp ; do i=1,id if (abs(tr_inp(i,j)-missing_value) > abs(roundoff*missing_value)) then - mask_in(i,j)=1.0 + mask_in(i,j) = 1.0 tr_inp(i,j) = tr_inp(i,j) * conversion else tr_inp(i,j) = missing_value @@ -837,7 +840,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t endif if ((G%mask2dT(i,j) == 1.0) .and. (z_edges_in(k) <= G%bathyT(i,j)) .and. & (mask_out(i,j) < 1.0)) & - fill(i,j)=1.0 + fill(i,j) = 1.0 enddo ; enddo call pass_var(fill, G%Domain) call pass_var(good, G%Domain) @@ -881,15 +884,15 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t enddo ! kd else - call time_interp_extern(fms_id, Time, data_in, verbose=.true., turns=turns) - do k=1,kd - do j=js,je - do i=is,ie - tr_z(i,j,k)=data_in(i,j,k) - if (abs(tr_z(i,j,k)-missing_value) < abs(roundoff*missing_value)) mask_z(i,j,k) = 0. - enddo + call time_interp_extern(fms_id, Time, data_in, verbose=.true., turns=turns) + do k=1,kd + do j=js,je + do i=is,ie + tr_z(i,j,k) = data_in(i,j,k) + if (abs(tr_z(i,j,k)-missing_value) < abs(roundoff*missing_value)) mask_z(i,j,k) = 0. enddo enddo + enddo endif end subroutine horiz_interp_and_extrap_tracer_fms_id @@ -901,9 +904,9 @@ subroutine meshgrid(x, y, x_T, y_T) real, dimension(size(x,1),size(y,1)), intent(inout) :: x_T !< output 2-dimensional array real, dimension(size(x,1),size(y,1)), intent(inout) :: y_T !< output 2-dimensional array - integer :: ni,nj,i,j + integer :: ni, nj, i, j - ni=size(x,1) ; nj=size(y,1) + ni = size(x,1) ; nj = size(y,1) do j=1,nj ; do i=1,ni x_T(i,j) = x(i) @@ -915,122 +918,4 @@ subroutine meshgrid(x, y, x_T, y_T) end subroutine meshgrid -! None of the subsequent code appears to be used at all. - -! !> Fill grid edges for integer data -! function fill_boundaries_int(m,cyclic_x,tripolar_n) result(mp) -! integer, dimension(:,:), intent(in) :: m !< input array (ND) -! logical, intent(in) :: cyclic_x !< True if domain is zonally re-entrant -! logical, intent(in) :: tripolar_n !< True if domain has an Arctic fold -! integer, dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp - -! real, dimension(size(m,1),size(m,2)) :: m_real -! real, dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp_real - -! m_real = real(m) - -! mp_real = fill_boundaries_real(m_real,cyclic_x,tripolar_n) - -! mp = int(mp_real) - -! end function fill_boundaries_int - -!> Fill grid edges for real data -! function fill_boundaries_real(m,cyclic_x,tripolar_n) result(mp) -! real, dimension(:,:), intent(in) :: m !< input array (ND) -! logical, intent(in) :: cyclic_x !< True if domain is zonally re-entrant -! logical, intent(in) :: tripolar_n !< True if domain has an Arctic fold -! real, dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp - -! integer :: ni,nj,i,j - -! ni=size(m,1); nj=size(m,2) - -! mp(1:ni,1:nj)=m(:,:) - -! if (cyclic_x) then -! mp(0,1:nj)=m(ni,1:nj) -! mp(ni+1,1:nj)=m(1,1:nj) -! else -! mp(0,1:nj)=m(1,1:nj) -! mp(ni+1,1:nj)=m(ni,1:nj) -! endif - -! mp(1:ni,0)=m(1:ni,1) -! if (tripolar_n) then -! do i=1,ni -! mp(i,nj+1)=m(ni-i+1,nj) -! enddo -! else -! mp(1:ni,nj+1)=m(1:ni,nj) -! endif - -! end function fill_boundaries_real - -!> Solve del2 (zi) = 0 using successive iterations -!! with a 5 point stencil. Only points fill==1 are -!! modified. Except where bad==1, information propagates -!! isotropically in index space. The resulting solution -!! in each region is an approximation to del2(zi)=0 subject to -!! boundary conditions along the valid points curve bounding this region. -! subroutine smooth_heights(zi,fill,bad,sor,niter,cyclic_x, tripolar_n) -! real, dimension(:,:), intent(inout) :: zi !< input and output array (ND) -! integer, dimension(size(zi,1),size(zi,2)), intent(in) :: fill !< same shape as zi, 1=fill -! integer, dimension(size(zi,1),size(zi,2)), intent(in) :: bad !< same shape as zi, 1=bad data -! real, intent(in) :: sor !< relaxation coefficient (ND) -! integer, intent(in) :: niter !< maximum number of iterations -! logical, intent(in) :: cyclic_x !< true if domain is zonally reentrant -! logical, intent(in) :: tripolar_n !< true if domain has an Arctic fold - -! ! Local variables -! real, dimension(size(zi,1),size(zi,2)) :: res, m -! integer, dimension(size(zi,1),size(zi,2),4) :: B -! real, dimension(0:size(zi,1)+1,0:size(zi,2)+1) :: mp -! integer, dimension(0:size(zi,1)+1,0:size(zi,2)+1) :: nm -! integer :: i,j,k,n -! integer :: ni,nj -! real :: Isum, bsum - -! ni=size(zi,1) ; nj=size(zi,2) - - -! mp(:,:) = fill_boundaries(zi,cyclic_x,tripolar_n) - -! B(:,:,:) = 0.0 -! nm(:,:) = fill_boundaries(bad,cyclic_x,tripolar_n) - -! do j=1,nj -! do i=1,ni -! if (fill(i,j) == 1) then -! B(i,j,1)=1-nm(i+1,j);B(i,j,2)=1-nm(i-1,j) -! B(i,j,3)=1-nm(i,j+1);B(i,j,4)=1-nm(i,j-1) -! endif -! enddo -! enddo - -! do n=1,niter -! do j=1,nj -! do i=1,ni -! if (fill(i,j) == 1) then -! bsum = real(B(i,j,1)+B(i,j,2)+B(i,j,3)+B(i,j,4)) -! Isum = 1.0/bsum -! res(i,j)=Isum*(B(i,j,1)*mp(i+1,j)+B(i,j,2)*mp(i-1,j)+& -! B(i,j,3)*mp(i,j+1)+B(i,j,4)*mp(i,j-1)) - mp(i,j) -! endif -! enddo -! enddo -! res(:,:)=res(:,:)*sor - -! do j=1,nj -! do i=1,ni -! mp(i,j)=mp(i,j)+res(i,j) -! enddo -! enddo - -! zi(:,:)=mp(1:ni,1:nj) -! mp = fill_boundaries(zi,cyclic_x,tripolar_n) -! enddo - -! end subroutine smooth_heights - end module MOM_horizontal_regridding From e4d984a7ad66527c3508e28c8f64798b76740044 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 14 Jan 2021 12:47:53 -0500 Subject: [PATCH 07/11] +Created MOM_diag_manager to wrap diag_manager Moved MOM_diag_manager_wrapper.F90 to MOM_diag_manager.F90 (mostly for brevity) and added use statements and provided interfaces for all of the diag_manager_mod, diag_data_mod and diag_axis_mod routines that are used within the MOM6 code, with some renaming of interfaces to reflect their use within the MOM6 code. All answers are bitwise identical. --- src/diagnostics/MOM_obsolete_diagnostics.F90 | 4 +- ...nager_wrapper.F90 => MOM_diag_manager.F90} | 19 +++++++--- src/framework/MOM_diag_mediator.F90 | 38 +++++++------------ src/framework/MOM_diag_remap.F90 | 4 +- src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 | 18 ++++----- 5 files changed, 39 insertions(+), 44 deletions(-) rename src/framework/{MOM_diag_manager_wrapper.F90 => MOM_diag_manager.F90} (88%) diff --git a/src/diagnostics/MOM_obsolete_diagnostics.F90 b/src/diagnostics/MOM_obsolete_diagnostics.F90 index e30749984d..bba8379bbb 100644 --- a/src/diagnostics/MOM_obsolete_diagnostics.F90 +++ b/src/diagnostics/MOM_obsolete_diagnostics.F90 @@ -4,10 +4,10 @@ module MOM_obsolete_diagnostics ! This file is part of MOM6. See LICENSE.md for the license. +use MOM_diag_manager, only : register_static_field_fms +use MOM_diag_mediator, only : diag_ctrl use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe use MOM_file_parser, only : param_file_type, log_version, get_param -use MOM_diag_mediator, only : diag_ctrl -use diag_manager_mod, only : register_static_field_fms=>register_static_field implicit none ; private diff --git a/src/framework/MOM_diag_manager_wrapper.F90 b/src/framework/MOM_diag_manager.F90 similarity index 88% rename from src/framework/MOM_diag_manager_wrapper.F90 rename to src/framework/MOM_diag_manager.F90 index 47dc701798..0c9f875bcd 100644 --- a/src/framework/MOM_diag_manager_wrapper.F90 +++ b/src/framework/MOM_diag_manager.F90 @@ -1,14 +1,23 @@ -!> A simple (very thin) wrapper for register_diag_field to avoid a compiler bug with PGI -module MOM_diag_manager_wrapper +!> A simple (very thin) wrapper for the FMS diag_manager routines, with some name changes +module MOM_diag_manager ! This file is part of MOM6. See LICENSE.md for the license. use MOM_time_manager, only : time_type +use diag_axis_mod, only : diag_axis_init, get_diag_axis_name, EAST, NORTH +use diag_data_mod, only : null_axis_id +use diag_manager_mod, only : diag_manager_init, diag_manager_end +use diag_manager_mod, only : send_data, diag_field_add_attribute, DIAG_FIELD_NOT_FOUND use diag_manager_mod, only : register_diag_field +use diag_manager_mod, only : register_static_field_fms=>register_static_field +use diag_manager_mod, only : get_diag_field_id_fms=>get_diag_field_id implicit none ; private -public register_diag_field_fms +public diag_manager_init, diag_manager_end +public diag_axis_init, get_diag_axis_name, EAST, NORTH, null_axis_id +public send_data, diag_field_add_attribute, DIAG_FIELD_NOT_FOUND +public register_diag_field_fms, register_static_field_fms, get_diag_field_id_fms !> A wrapper for register_diag_field_array() interface register_diag_field_fms @@ -85,11 +94,11 @@ integer function register_diag_field_scalar_fms(module_name, field_name, init_ti end function register_diag_field_scalar_fms -!> \namespace mom_diag_manager_wrapper +!> \namespace mom_diag_manager !! !! This module simply wraps register_diag_field() from FMS's diag_manager_mod. !! We used to be able to import register_diag_field and rename it to register_diag_field_fms !! with a simple "use, only : register_diag_field_fms => register_diag_field" but PGI 16.5 !! has a bug that refuses to compile this - earlier versions did work. -end module MOM_diag_manager_wrapper +end module MOM_diag_manager diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index fa4a4a2701..071585a951 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -9,37 +9,28 @@ module MOM_diag_mediator use MOM_coms, only : PE_here use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_MODULE, CLOCK_ROUTINE +use MOM_diag_manager, only : diag_manager_init, diag_manager_end +use MOM_diag_manager, only : diag_axis_init, get_diag_axis_name, null_axis_id +use MOM_diag_manager, only : send_data, diag_field_add_attribute, EAST, NORTH +use MOM_diag_manager, only : register_diag_field_fms, register_static_field_fms +use MOM_diag_manager, only : get_diag_field_id_fms, DIAG_FIELD_NOT_FOUND +use MOM_diag_remap, only : diag_remap_ctrl, diag_remap_update, diag_remap_calc_hmask +use MOM_diag_remap, only : diag_remap_init, diag_remap_end, diag_remap_do_remap +use MOM_diag_remap, only : vertically_reintegrate_diag_field, vertically_interpolate_diag_field +use MOM_diag_remap, only : horizontally_average_diag_field, diag_remap_get_axes_info +use MOM_diag_remap, only : diag_remap_configure_axes, diag_remap_axes_configured +use MOM_diag_remap, only : diag_remap_diag_registration_closed, diag_remap_set_active +use MOM_EOS, only : EOS_type use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe, assert use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type -use MOM_io, only : slasher, vardesc, query_vardesc, mom_read_data +use MOM_io, only : slasher, vardesc, query_vardesc, MOM_read_data use MOM_io, only : get_filename_appendix use MOM_safe_alloc, only : safe_alloc_ptr, safe_alloc_alloc use MOM_string_functions, only : lowercase use MOM_time_manager, only : time_type use MOM_unit_scaling, only : unit_scale_type use MOM_verticalGrid, only : verticalGrid_type -use MOM_EOS, only : EOS_type -use MOM_diag_remap, only : diag_remap_ctrl -use MOM_diag_remap, only : diag_remap_update -use MOM_diag_remap, only : diag_remap_calc_hmask -use MOM_diag_remap, only : diag_remap_init, diag_remap_end, diag_remap_do_remap -use MOM_diag_remap, only : vertically_reintegrate_diag_field, vertically_interpolate_diag_field -use MOM_diag_remap, only : diag_remap_configure_axes, diag_remap_axes_configured -use MOM_diag_remap, only : diag_remap_get_axes_info, diag_remap_set_active -use MOM_diag_remap, only : diag_remap_diag_registration_closed -use MOM_diag_remap, only : horizontally_average_diag_field - -use diag_axis_mod, only : get_diag_axis_name -use diag_data_mod, only : null_axis_id -use diag_manager_mod, only : diag_manager_init, diag_manager_end -use diag_manager_mod, only : send_data, diag_axis_init, EAST, NORTH, diag_field_add_attribute -! The following module is needed for PGI since the following line does not compile with PGI 6.5.0 -! was: use diag_manager_mod, only : register_diag_field_fms=>register_diag_field -use MOM_diag_manager_wrapper, only : register_diag_field_fms -use diag_manager_mod, only : register_static_field_fms=>register_static_field -use diag_manager_mod, only : get_diag_field_id_fms=>get_diag_field_id -use diag_manager_mod, only : DIAG_FIELD_NOT_FOUND implicit none ; private @@ -482,10 +473,9 @@ subroutine set_axes_info(G, GV, US, param_file, diag_cs, set_vertical) call define_axes_group(diag_cs, (/ id_xh, id_yq /), diag_cs%axesCv1, & x_cell_method='mean', y_cell_method='point', is_v_point=.true.) - ! Axis group for special null axis from diag manager + ! Axis group for special null axis from diag manager. (Could null_axis_id be made MOM specific?) call define_axes_group(diag_cs, (/ null_axis_id /), diag_cs%axesNull) - !Non-native Non-downsampled if (diag_cs%num_diag_coords>0) then allocate(diag_cs%remap_axesZL(diag_cs%num_diag_coords)) diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index f27a153a2b..462ba7bf5e 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -60,6 +60,8 @@ module MOM_diag_remap use MOM_coms, only : reproducing_sum_EFP, EFP_to_real use MOM_coms, only : EFP_type, assignment(=), EFP_sum_across_PEs use MOM_error_handler, only : MOM_error, FATAL, assert, WARNING +use MOM_debugging, only : check_column_integrals +use MOM_diag_manager, only : diag_axis_init use MOM_diag_vkernels, only : interpolate_column, reintegrate_column use MOM_file_parser, only : get_param, log_param, param_file_type use MOM_io, only : slasher, mom_read_data @@ -80,9 +82,7 @@ module MOM_diag_remap use coord_sigma, only : build_sigma_column use coord_rho, only : build_rho_column -use diag_manager_mod, only : diag_axis_init -use MOM_debugging, only : check_column_integrals implicit none ; private public diag_remap_ctrl diff --git a/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 b/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 index 90ae47450d..397696c0ba 100644 --- a/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 +++ b/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 @@ -3,19 +3,15 @@ module MOM_IS_diag_mediator ! This file is a part of SIS2. See LICENSE.md for the license. -use MOM_grid, only : ocean_grid_type - -use MOM_coms, only : PE_here +use MOM_coms, only : PE_here +use MOM_diag_manager, only : diag_manager_init, send_data, diag_axis_init, EAST, NORTH +use MOM_diag_manager, only : register_diag_field_fms, register_static_field_fms use MOM_error_handler, only : MOM_error, FATAL, is_root_pe, assert -use MOM_file_parser, only : get_param, log_param, log_version, param_file_type -use MOM_safe_alloc, only : safe_alloc_ptr, safe_alloc_alloc +use MOM_file_parser, only : get_param, log_param, log_version, param_file_type +use MOM_grid, only : ocean_grid_type +use MOM_safe_alloc, only : safe_alloc_ptr, safe_alloc_alloc use MOM_string_functions, only : lowercase, uppercase, slasher -use MOM_time_manager, only : time_type - -use diag_manager_mod, only : diag_manager_init -use diag_manager_mod, only : send_data, diag_axis_init,EAST,NORTH -use diag_manager_mod, only : register_diag_field_fms=>register_diag_field -use diag_manager_mod, only : register_static_field_fms=>register_static_field +use MOM_time_manager, only : time_type implicit none ; private From 76b3ccce48a1d4fb1c61fe6aee0a06d5b60a6219 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 14 Jan 2021 13:57:13 -0500 Subject: [PATCH 08/11] Use only for netcdf in MOM_horizontal_regridding Made netcdf module use dependencies explicit in MOM_horizontal_regridding.F90. All answers are bitwise identical. --- src/framework/MOM_horizontal_regridding.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index f1ab073938..b91509cc1d 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -16,7 +16,8 @@ module MOM_horizontal_regridding use MOM_io_wrapper, only : axistype, get_axis_data use MOM_time_manager, only : time_type -use netcdf +use netcdf, only : NF90_OPEN, NF90_NOWRITE, NF90_GET_ATT, NF90_GET_VAR +use netcdf, only : NF90_INQ_VARID, NF90_INQUIRE_VARIABLE, NF90_INQUIRE_DIMENSION implicit none ; private From e29a12b80a5fbcded3ee50ffa26f2ec1dd0d9ff2 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Fri, 15 Jan 2021 15:19:19 -0500 Subject: [PATCH 09/11] Adds a GH workflow to check driver APIs - Checks that MOM6 can be compiled with the various coupled drivers. --- .github/workflows/coupled-api.yml | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) create mode 100644 .github/workflows/coupled-api.yml diff --git a/.github/workflows/coupled-api.yml b/.github/workflows/coupled-api.yml new file mode 100644 index 0000000000..4380f102ea --- /dev/null +++ b/.github/workflows/coupled-api.yml @@ -0,0 +1,30 @@ +name: API for coupled drivers + +on: [push, pull_request] + +jobs: + test-top-api: + + runs-on: ubuntu-latest + defaults: + run: + working-directory: .testing + + steps: + - uses: actions/checkout@v2 + with: + submodules: recursive + + - uses: ./.github/actions/testing-setup + + - name: Compile MOM6 for the GFDL coupled driver + shell: bash + run: make check_mom6_api_coupled -j + + - name: Compile MOM6 for the NUOPC driver + shell: bash + run: make check_mom6_api_nuopc -j + + - name: Compile MOM6 for the MCT driver + shell: bash + run: make check_mom6_api_mct -j From eb0c03a3dd80c04f69a6b685ffb1553f20ed04e0 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Fri, 15 Jan 2021 16:09:55 -0500 Subject: [PATCH 10/11] Avoid unnecessary steps in actions/testing_setup - We spent 1 minute build symmetric unnecessarily when checking the coupled APIs --- .github/actions/testing-setup/action.yml | 13 +++++++++++-- .github/workflows/coupled-api.yml | 3 +++ 2 files changed, 14 insertions(+), 2 deletions(-) diff --git a/.github/actions/testing-setup/action.yml b/.github/actions/testing-setup/action.yml index c6fae4ad58..c47270af0d 100644 --- a/.github/actions/testing-setup/action.yml +++ b/.github/actions/testing-setup/action.yml @@ -1,5 +1,14 @@ name: 'Build-.testing-prerequisites' description: 'Build pre-requisites for .testing including FMS and a symmetric MOM6 executable' +inputs: + build_symmetric: + description: 'If true, will build the symmetric MOM6 executable' + required: false + default: 'true' + install_python: + description: 'If true, will install the local python env needed for .testing' + required: false + default: 'true' runs: using: 'composite' steps: @@ -51,7 +60,7 @@ runs: run: | echo "::group::Compile MOM6 in symmetric memory mode" cd .testing - make build/symmetric/MOM6 -j + test ${{ inputs.build_symmetric }} == true && make build/symmetric/MOM6 -j echo "::endgroup::" - name: Install local python venv for generating input data @@ -59,7 +68,7 @@ runs: run: | echo "::group::Create local python env for input data generation" cd .testing - make work/local-env + test ${{ inputs.install_python }} == true && make work/local-env echo "::endgroup::" - name: Set flags diff --git a/.github/workflows/coupled-api.yml b/.github/workflows/coupled-api.yml index 4380f102ea..86d7262548 100644 --- a/.github/workflows/coupled-api.yml +++ b/.github/workflows/coupled-api.yml @@ -16,6 +16,9 @@ jobs: submodules: recursive - uses: ./.github/actions/testing-setup + with: + build_symmetric: 'false' + install_python: 'false' - name: Compile MOM6 for the GFDL coupled driver shell: bash From 354d6d6406fc09f534e954c6ba42881a6116e725 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Mon, 18 Jan 2021 09:00:59 -0500 Subject: [PATCH 11/11] MOM_hor_visc: Revert tension/shear loop fusion The fusion of tension and shear strains yields a 1-2% speedup, but also breaks the style convention of capitalization of vertex points, and also evaluates the tension terms over a slightly larger domain, so it has been reverted. A note has been added to investigate this later. --- .../lateral/MOM_hor_visc.F90 | 25 ++++++++++++------- 1 file changed, 16 insertions(+), 9 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 06a6269dc5..003b134b2a 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -528,21 +528,29 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! shearing strain advocated by Smagorinsky (1993) and discussed in ! Griffies and Hallberg (2000). - do j=Jsq-2,Jeq+2 ; do i=Isq-2,Ieq+2 - ! Calculate horizontal tension + ! NOTE: There is a ~1% speedup when the tension and shearing loops below + ! are fused (presumably due to shared access of Id[xy]C[uv]). However, + ! this breaks the center/vertex index case convention, and also evaluates + ! the dudx and dvdy terms beyond their valid bounds. + ! TODO: Explore methods for retaining both the syntax and speedup. + + ! Calculate horizontal tension + do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2 dudx(i,j) = CS%DY_dxT(i,j)*(G%IdyCu(I,j) * u(I,j,k) - & G%IdyCu(I-1,j) * u(I-1,j,k)) dvdy(i,j) = CS%DX_dyT(i,j)*(G%IdxCv(i,J) * v(i,J,k) - & G%IdxCv(i,J-1) * v(i,J-1,k)) sh_xx(i,j) = dudx(i,j) - dvdy(i,j) + enddo ; enddo - ! Components for the shearing strain + ! Components for the shearing strain + do J=Jsq-2,Jeq+2 ; do I=Isq-2,Ieq+2 dvdx(I,J) = CS%DY_dxBu(I,J)*(v(i+1,J,k)*G%IdyCv(i+1,J) - v(i,J,k)*G%IdyCv(i,J)) dudy(I,J) = CS%DX_dyBu(I,J)*(u(I,j+1,k)*G%IdxCu(I,j+1) - u(I,j,k)*G%IdxCu(I,j)) enddo ; enddo if (CS%id_normstress > 0) then - do j=Jsq-2,Jeq+2 ; do i=Isq-2,Ieq+2 + do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2 NoSt(i,j,k) = sh_xx(i,j) enddo ; enddo endif @@ -885,6 +893,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, Kh(i,j) = CS%Kh_bg_xx(i,j) enddo ; enddo + ! NOTE: The following do-block can be decomposed and vectorized after the + ! stack size has been reduced. do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 if (CS%add_LES_viscosity) then if (CS%Smagorinsky_Kh) & @@ -1217,12 +1227,9 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! All viscosity contributions above are subject to resolution scaling + ! NOTE: The following do-block can be decomposed and vectorized after the + ! stack size has been reduced. do J=js-1,Jeq ; do I=is-1,Ieq - ! NOTE: The following do-block can be decomposed and vectorized, but - ! appears to cause slowdown on some machines. Evidence suggests that - ! this is caused by excessive spilling of stack variables. - ! TODO: Vectorize these loops after stack usage has been reduced.. - if (rescale_Kh) & Kh(i,j) = VarMix%Res_fn_q(i,j) * Kh(i,j)