diff --git a/.travis.yml b/.travis.yml index efcd7344c9..4d608c9606 100644 --- a/.travis.yml +++ b/.travis.yml @@ -54,11 +54,12 @@ script: - (cd MOM6-examples/ocean_only/benchmark && mkdir -p RESTART && mpirun -n 1 ../../../build/ocn/MOM6 1> log) - (cd MOM6-examples/ocean_only/circle_obcs && mkdir -p RESTART && mpirun -n 1 ../../../build/symocn/MOM6 1> log) - md5sum MOM6-examples/ocean_only/{*,*/*}/ocean.stats > stats.md5 + - (cd MOM6-examples/ocean_only/unit_tests && mkdir -p RESTART && mpirun -n 1 ../../../build/ocn/MOM6 1> log) - (cd MOM6-examples/ocean_only/double_gyre && mkdir -p RESTART && mpirun -n 2 ../../../build/symocn/MOM6 1> log) - - (cd MOM6-examples/ocean_only/flow_downslope/layer && mkdir -p RESTART && mpirun -n 2 ../../../../build/ocn/MOM6 1>log) - - (cd MOM6-examples/ocean_only/flow_downslope/rho && mkdir -p RESTART && mpirun -n 2 ../../../../build/ocn/MOM6 1>log) - - (cd MOM6-examples/ocean_only/flow_downslope/sigma && mkdir -p RESTART && mpirun -n 2 ../../../../build/ocn/MOM6 1>log) - - (cd MOM6-examples/ocean_only/flow_downslope/z && mkdir -p RESTART && mpirun -n 2 ../../../../build/ocn/MOM6 1>log) + - (cd MOM6-examples/ocean_only/flow_downslope/layer && mkdir -p RESTART && mpirun -n 2 ../../../../build/symocn/MOM6 1>log) + - (cd MOM6-examples/ocean_only/flow_downslope/rho && mkdir -p RESTART && mpirun -n 2 ../../../../build/symocn/MOM6 1>log) + - (cd MOM6-examples/ocean_only/flow_downslope/sigma && mkdir -p RESTART && mpirun -n 2 ../../../../build/symocn/MOM6 1>log) + - (cd MOM6-examples/ocean_only/flow_downslope/z && mkdir -p RESTART && mpirun -n 2 ../../../../build/symocn/MOM6 1>log) - (cd MOM6-examples/ocean_only/benchmark && mkdir -p RESTART && mpirun -n 4 ../../../build/symocn/MOM6 1> log) - (cd MOM6-examples/ocean_only/circle_obcs && mkdir -p RESTART && mpirun -n 4 ../../../build/symocn/MOM6 1> log) - md5sum -c stats.md5 diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index 62a85148e8..75624104c3 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -95,7 +95,7 @@ module MOM_ALE integer, dimension(:), allocatable :: id_Htracer_remap_tendency !< diagnostic id integer, dimension(:), allocatable :: id_Htracer_remap_tendency_2d !< diagnostic id logical, dimension(:), allocatable :: do_tendency_diag !< flag for doing diagnostics - integer :: id_dzRegrid + integer :: id_dzRegrid = -1 !< diagnostic id end type @@ -669,7 +669,7 @@ subroutine remap_all_state_vars(CS_remapping, CS_ALE, G, GV, h_old, h_new, dxInt ! Build the start and final grids h1(:) = h_old(i,j,:) h2(:) = h_new(i,j,:) - call remapping_core_h( nz, h1, Reg%Tr(m)%t(i,j,:), nz, h2, u_column, CS_remapping ) + call remapping_core_h(CS_remapping, nz, h1, Reg%Tr(m)%t(i,j,:), nz, h2, u_column) ! Intermediate steps for tendency of tracer concentration and tracer content. ! Note: do not merge the two if-tests, since do_tendency_diag(:) is not @@ -764,7 +764,7 @@ subroutine remap_all_state_vars(CS_remapping, CS_ALE, G, GV, h_old, h_new, dxInt else h2(:) = 0.5 * ( h_new(i,j,:) + h_new(i+1,j,:) ) endif - call remapping_core_h( nz, h1, u(I,j,:), nz, h2, u_column, CS_remapping ) + call remapping_core_h(CS_remapping, nz, h1, u(I,j,:), nz, h2, u_column) u(I,j,:) = u_column(:) endif enddo @@ -789,7 +789,7 @@ subroutine remap_all_state_vars(CS_remapping, CS_ALE, G, GV, h_old, h_new, dxInt else h2(:) = 0.5 * ( h_new(i,j,:) + h_new(i,j+1,:) ) endif - call remapping_core_h( nz, h1, v(i,J,:), nz, h2, u_column, CS_remapping ) + call remapping_core_h(CS_remapping, nz, h1, v(i,J,:), nz, h2, u_column) v(i,J,:) = u_column(:) endif enddo @@ -849,7 +849,7 @@ subroutine ALE_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_c call dzFromH1H2( n_points, h_src(i,j,1:n_points), GV%ke, h_dst(i,j,:), dx ) call remapping_core_w(CS, n_points, h_src(i,j,1:n_points), s_src(i,j,1:n_points), GV%ke, dx, s_dst(i,j,:)) else - call remapping_core_h(n_points, h_src(i,j,1:n_points), s_src(i,j,1:n_points), GV%ke, h_dst(i,j,:), s_dst(i,j,:), CS) + call remapping_core_h(CS, n_points, h_src(i,j,1:n_points), s_src(i,j,1:n_points), GV%ke, h_dst(i,j,:), s_dst(i,j,:)) endif else s_dst(i,j,:) = 0. diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index e9726ae06a..b6909e7fa2 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -76,7 +76,7 @@ module MOM_remapping !! changing the numerical result, except where !! a division by zero would otherwise occur. -logical, parameter :: old_algorithm = .true. !< Use the old "broken" algorithm. +logical, parameter :: old_algorithm = .false. !< Use the old "broken" algorithm. !! This is a temporary measure to assist !! debugging until we delete the old algorithm. @@ -159,14 +159,14 @@ end function isPosSumErrSignificant !> Remaps column of values u0 on grid h0 to grid h1 !! assuming the top edge is aligned. -subroutine remapping_core_h( n0, h0, u0, n1, h1, u1, CS ) +subroutine remapping_core_h(CS, n0, h0, u0, n1, h1, u1) + type(remapping_CS), intent(in) :: CS !< Remapping control structure integer, intent(in) :: n0 !< Number of cells on source grid real, dimension(n0), intent(in) :: h0 !< Cell widths on source grid real, dimension(n0), intent(in) :: u0 !< Cell averages on source grid integer, intent(in) :: n1 !< Number of cells on target grid real, dimension(n1), intent(in) :: h1 !< Cell widths on target grid real, dimension(n1), intent(out) :: u1 !< Cell averages on target grid - type(remapping_CS), intent(in) :: CS !< Remapping control structure ! Local variables integer :: iMethod real, dimension(n0,2) :: ppoly_r_E !Edge value of polynomial @@ -628,6 +628,8 @@ subroutine remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefficients, n1, h isrc_max(i0) = i_max i_max = i_sub + 1 dh_max = 0. + ! Record the source cell thickness found by summing the sub-cell thicknesses. + h0_eff(i0) = dh0_eff if (i0 < n0) then i0 = i0 + 1 h0_supply = h0(i0) @@ -1524,7 +1526,9 @@ end subroutine end_remapping !> Runs unit tests on remapping functions. !! Should only be called from a single/root thread !! Returns True if a test fails, otherwise False -logical function remapping_unit_tests() +logical function remapping_unit_tests(verbose) + logical, intent(in) :: verbose !< If true, write results to stdout + ! Local variables integer, parameter :: n0 = 4, n1 = 3, n2 = 6 real :: h0(n0), x0(n0+1), u0(n0) real :: h1(n1), x1(n1+1), u1(n1), hn1(n1), dx1(n1+1) @@ -1537,9 +1541,11 @@ logical function remapping_unit_tests() real, allocatable, dimension(:,:) :: ppoly0_E, ppoly0_S, ppoly0_coefficients integer :: i real :: err - logical :: thisTest + logical :: thisTest, v - write(*,*) '===== MOM_remapping: remapping_unit_tests ================' + v = verbose + + write(*,*) '==== MOM_remapping: remapping_unit_tests =================' remapping_unit_tests = .false. ! Normally return false thisTest = .false. @@ -1560,8 +1566,8 @@ logical function remapping_unit_tests() thisTest = .false. call initialize_remapping(CS, 'PPM_H4') - write(*,*) 'h0 (test data)' - call dumpGrid(n0,h0,x0,u0) + if (verbose) write(*,*) 'h0 (test data)' + if (verbose) call dumpGrid(n0,h0,x0,u0) call dzFromH1H2( n0, h0, n1, h1, dx1 ) call remapping_core_w( CS, n0, h0, u0, n1, dx1, u1 ) @@ -1569,8 +1575,8 @@ logical function remapping_unit_tests() err=u1(i)-8.*(0.5*real(1+n1)-real(i)) if (abs(err)>real(n1-1)*epsilon(err)) thisTest = .true. enddo - write(*,*) 'h1 (by projection)' - call dumpGrid(n1,h1,x1,u1) + if (verbose) write(*,*) 'h1 (by projection)' + if (verbose) call dumpGrid(n1,h1,x1,u1) if (thisTest) write(*,*) 'remapping_unit_tests: Failed remapping_core_w()' remapping_unit_tests = remapping_unit_tests .or. thisTest @@ -1601,8 +1607,8 @@ logical function remapping_unit_tests() call remapByDeltaZ( n0, h0, u0, ppoly0_E, ppoly0_coefficients, & n1, x1-x0(1:n1+1), & INTEGRATION_PPM, u1, hn1 ) - write(*,*) 'h1 (by delta)' - call dumpGrid(n1,h1,x1,u1) + if (verbose) write(*,*) 'h1 (by delta)' + if (verbose) call dumpGrid(n1,h1,x1,u1) hn1=hn1-h1 do i=1,n1 err=u1(i)-8.*(0.5*real(1+n1)-real(i)) @@ -1618,10 +1624,10 @@ logical function remapping_unit_tests() call remapByDeltaZ( n0, h0, u0, ppoly0_E, ppoly0_coefficients, & n2, dx2, & INTEGRATION_PPM, u2, hn2 ) - write(*,*) 'h2' - call dumpGrid(n2,h2,x2,u2) - write(*,*) 'hn2' - call dumpGrid(n2,hn2,x2,u2) + if (verbose) write(*,*) 'h2' + if (verbose) call dumpGrid(n2,h2,x2,u2) + if (verbose) write(*,*) 'hn2' + if (verbose) call dumpGrid(n2,hn2,x2,u2) do i=1,n2 err=u2(i)-8./2.*(0.5*real(1+n2)-real(i)) @@ -1630,11 +1636,11 @@ logical function remapping_unit_tests() if (thisTest) write(*,*) 'remapping_unit_tests: Failed remapByDeltaZ() 2' remapping_unit_tests = remapping_unit_tests .or. thisTest - write(*,*) 'Via sub-cells' + if (verbose) write(*,*) 'Via sub-cells' thisTest = .false. call remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefficients, & n2, h2, INTEGRATION_PPM, .false., u2, err ) - call dumpGrid(n2,h2,x2,u2) + if (verbose) call dumpGrid(n2,h2,x2,u2) do i=1,n2 err=u2(i)-8./2.*(0.5*real(1+n2)-real(i)) @@ -1645,12 +1651,13 @@ logical function remapping_unit_tests() call remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefficients, & 6, (/.125,.125,.125,.125,.125,.125/), INTEGRATION_PPM, .false., u2, err ) - call dumpGrid(6,h2,x2,u2) + if (verbose) call dumpGrid(6,h2,x2,u2) call remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefficients, & 3, (/2.25,1.5,1./), INTEGRATION_PPM, .false., u2, err ) - call dumpGrid(3,h2,x2,u2) + if (verbose) call dumpGrid(3,h2,x2,u2) + if (.not. remapping_unit_tests) write(*,*) 'Pass' write(*,*) '===== MOM_remapping: new remapping_unit_tests ==================' @@ -1660,101 +1667,107 @@ logical function remapping_unit_tests() allocate(ppoly0_S(5,2)) call PCM_reconstruction(3, (/1.,2.,4./), ppoly0_E(1:3,:), ppoly0_coefficients(1:3,:) ) - remapping_unit_tests = remapping_unit_tests .or. test_answer(3, ppoly0_E(:,1), (/1.,2.,4./), 'PCM: left edges') - remapping_unit_tests = remapping_unit_tests .or. test_answer(3, ppoly0_E(:,2), (/1.,2.,4./), 'PCM: right edges') - remapping_unit_tests = remapping_unit_tests .or. test_answer(3, ppoly0_coefficients(:,1), (/1.,2.,4./), 'PCM: P0') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 3, ppoly0_E(:,1), (/1.,2.,4./), 'PCM: left edges') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 3, ppoly0_E(:,2), (/1.,2.,4./), 'PCM: right edges') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 3, ppoly0_coefficients(:,1), (/1.,2.,4./), 'PCM: P0') call PLM_reconstruction(3, (/1.,1.,1./), (/1.,3.,5./), ppoly0_E(1:3,:), ppoly0_coefficients(1:3,:) ) - remapping_unit_tests = remapping_unit_tests .or. test_answer(3, ppoly0_E(:,1), (/1.,2.,5./), 'Unlim PLM: left edges') - remapping_unit_tests = remapping_unit_tests .or. test_answer(3, ppoly0_E(:,2), (/1.,4.,5./), 'Unlim PLM: right edges') - remapping_unit_tests = remapping_unit_tests .or. test_answer(3, ppoly0_coefficients(:,1), (/1.,2.,5./), 'Unlim PLM: P0') - remapping_unit_tests = remapping_unit_tests .or. test_answer(3, ppoly0_coefficients(:,2), (/0.,2.,0./), 'Unlim PLM: P1') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 3, ppoly0_E(:,1), (/1.,2.,5./), 'Unlim PLM: left edges') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 3, ppoly0_E(:,2), (/1.,4.,5./), 'Unlim PLM: right edges') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 3, ppoly0_coefficients(:,1), (/1.,2.,5./), 'Unlim PLM: P0') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 3, ppoly0_coefficients(:,2), (/0.,2.,0./), 'Unlim PLM: P1') call PLM_reconstruction(3, (/1.,1.,1./), (/1.,2.,7./), ppoly0_E(1:3,:), ppoly0_coefficients(1:3,:) ) - remapping_unit_tests = remapping_unit_tests .or. test_answer(3, ppoly0_E(:,1), (/1.,1.,7./), 'Left lim PLM: left edges') - remapping_unit_tests = remapping_unit_tests .or. test_answer(3, ppoly0_E(:,2), (/1.,3.,7./), 'Left lim PLM: right edges') - remapping_unit_tests = remapping_unit_tests .or. test_answer(3, ppoly0_coefficients(:,1), (/1.,1.,7./), 'Left lim PLM: P0') - remapping_unit_tests = remapping_unit_tests .or. test_answer(3, ppoly0_coefficients(:,2), (/0.,2.,0./), 'Left lim PLM: P1') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 3, ppoly0_E(:,1), (/1.,1.,7./), 'Left lim PLM: left edges') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 3, ppoly0_E(:,2), (/1.,3.,7./), 'Left lim PLM: right edges') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 3, ppoly0_coefficients(:,1), (/1.,1.,7./), 'Left lim PLM: P0') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 3, ppoly0_coefficients(:,2), (/0.,2.,0./), 'Left lim PLM: P1') call PLM_reconstruction(3, (/1.,1.,1./), (/1.,6.,7./), ppoly0_E(1:3,:), ppoly0_coefficients(1:3,:) ) - remapping_unit_tests = remapping_unit_tests .or. test_answer(3, ppoly0_E(:,1), (/1.,5.,7./), 'Right lim PLM: left edges') - remapping_unit_tests = remapping_unit_tests .or. test_answer(3, ppoly0_E(:,2), (/1.,7.,7./), 'Right lim PLM: right edges') - remapping_unit_tests = remapping_unit_tests .or. test_answer(3, ppoly0_coefficients(:,1), (/1.,5.,7./), 'Right lim PLM: P0') - remapping_unit_tests = remapping_unit_tests .or. test_answer(3, ppoly0_coefficients(:,2), (/0.,2.,0./), 'Right lim PLM: P1') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 3, ppoly0_E(:,1), (/1.,5.,7./), 'Right lim PLM: left edges') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 3, ppoly0_E(:,2), (/1.,7.,7./), 'Right lim PLM: right edges') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 3, ppoly0_coefficients(:,1), (/1.,5.,7./), 'Right lim PLM: P0') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 3, ppoly0_coefficients(:,2), (/0.,2.,0./), 'Right lim PLM: P1') call PLM_reconstruction(3, (/1.,2.,3./), (/1.,4.,9./), ppoly0_E(1:3,:), ppoly0_coefficients(1:3,:) ) - remapping_unit_tests = remapping_unit_tests .or. test_answer(3, ppoly0_E(:,1), (/1.,2.,9./), 'Non-uniform line PLM: left edges') - remapping_unit_tests = remapping_unit_tests .or. test_answer(3, ppoly0_E(:,2), (/1.,6.,9./), 'Non-uniform line PLM: right edges') - remapping_unit_tests = remapping_unit_tests .or. test_answer(3, ppoly0_coefficients(:,1), (/1.,2.,9./), 'Non-uniform line PLM: P0') - remapping_unit_tests = remapping_unit_tests .or. test_answer(3, ppoly0_coefficients(:,2), (/0.,4.,0./), 'Non-uniform line PLM: P1') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 3, ppoly0_E(:,1), (/1.,2.,9./), 'Non-uniform line PLM: left edges') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 3, ppoly0_E(:,2), (/1.,6.,9./), 'Non-uniform line PLM: right edges') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 3, ppoly0_coefficients(:,1), (/1.,2.,9./), 'Non-uniform line PLM: P0') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 3, ppoly0_coefficients(:,2), (/0.,4.,0./), 'Non-uniform line PLM: P1') call edge_values_explicit_h4( 5, (/1.,1.,1.,1.,1./), (/1.,3.,5.,7.,9./), ppoly0_E ) - remapping_unit_tests = remapping_unit_tests .or. test_answer(5, ppoly0_E(:,1), (/0.,2.,4.,6.,8./), 'Line H4: left edges') - remapping_unit_tests = remapping_unit_tests .or. test_answer(5, ppoly0_E(:,2), (/2.,4.,6.,8.,10./), 'Line H4: right edges') - remapping_unit_tests = .false. ! edge_values_explicit_h4 fails due to roundoff! + thisTest = test_answer(v, 5, ppoly0_E(:,1), (/0.,2.,4.,6.,8./), 'Line H4: left edges') ! Currently fails due to roundoff + thisTest = test_answer(v, 5, ppoly0_E(:,2), (/2.,4.,6.,8.,10./), 'Line H4: right edges') ! Currently fails due to roundoff ppoly0_E(:,1) = (/0.,2.,4.,6.,8./) ppoly0_E(:,2) = (/2.,4.,6.,8.,10./) call PPM_reconstruction(5, (/1.,1.,1.,1.,1./), (/1.,3.,5.,7.,9./), ppoly0_E(1:5,:), ppoly0_coefficients(1:5,:) ) - remapping_unit_tests = remapping_unit_tests .or. test_answer(5, ppoly0_coefficients(:,1), (/1.,2.,4.,6.,9./), 'Line PPM: P0') - remapping_unit_tests = remapping_unit_tests .or. test_answer(5, ppoly0_coefficients(:,2), (/0.,2.,2.,2.,0./), 'Line PPM: P1') - remapping_unit_tests = remapping_unit_tests .or. test_answer(5, ppoly0_coefficients(:,3), (/0.,0.,0.,0.,0./), 'Line PPM: P2') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 5, ppoly0_coefficients(:,1), (/1.,2.,4.,6.,9./), 'Line PPM: P0') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 5, ppoly0_coefficients(:,2), (/0.,2.,2.,2.,0./), 'Line PPM: P1') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 5, ppoly0_coefficients(:,3), (/0.,0.,0.,0.,0./), 'Line PPM: P2') call edge_values_explicit_h4( 5, (/1.,1.,1.,1.,1./), (/1.,1.,7.,19.,37./), ppoly0_E ) - remapping_unit_tests = remapping_unit_tests .or. test_answer(5, ppoly0_E(:,1), (/3.,0.,3.,12.,27./), 'Parabola H4: left edges') - remapping_unit_tests = remapping_unit_tests .or. test_answer(5, ppoly0_E(:,2), (/0.,3.,12.,27.,48./), 'Parabola H4: right edges') - remapping_unit_tests = .false. ! edge_values_explicit_h4 fails due to roundoff! + thisTest = test_answer(v, 5, ppoly0_E(:,1), (/3.,0.,3.,12.,27./), 'Parabola H4: left edges') ! Currently fails due to roundoff + thisTest = test_answer(v, 5, ppoly0_E(:,2), (/0.,3.,12.,27.,48./), 'Parabola H4: right edges') ! Currently fails due to roundoff ppoly0_E(:,1) = (/0.,0.,3.,12.,27./) ppoly0_E(:,2) = (/0.,3.,12.,27.,48./) call PPM_reconstruction(5, (/1.,1.,1.,1.,1./), (/0.,1.,7.,19.,37./), ppoly0_E(1:5,:), ppoly0_coefficients(1:5,:) ) - remapping_unit_tests = remapping_unit_tests .or. test_answer(5, ppoly0_E(:,1), (/0.,0.,3.,12.,37./), 'Parabola PPM: left edges') - remapping_unit_tests = remapping_unit_tests .or. test_answer(5, ppoly0_E(:,2), (/0.,3.,12.,27.,37./), 'Parabola PPM: right edges') - remapping_unit_tests = remapping_unit_tests .or. test_answer(5, ppoly0_coefficients(:,1), (/0.,0.,3.,12.,37./), 'Parabola PPM: P0') - remapping_unit_tests = remapping_unit_tests .or. test_answer(5, ppoly0_coefficients(:,2), (/0.,0.,6.,12.,0./), 'Parabola PPM: P1') - remapping_unit_tests = remapping_unit_tests .or. test_answer(5, ppoly0_coefficients(:,3), (/0.,3.,3.,3.,0./), 'Parabola PPM: P2') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 5, ppoly0_E(:,1), (/0.,0.,3.,12.,37./), 'Parabola PPM: left edges') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 5, ppoly0_E(:,2), (/0.,3.,12.,27.,37./), 'Parabola PPM: right edges') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 5, ppoly0_coefficients(:,1), (/0.,0.,3.,12.,37./), 'Parabola PPM: P0') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 5, ppoly0_coefficients(:,2), (/0.,0.,6.,12.,0./), 'Parabola PPM: P1') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 5, ppoly0_coefficients(:,3), (/0.,3.,3.,3.,0./), 'Parabola PPM: P2') ppoly0_E(:,1) = (/0.,0.,6.,10.,15./) ppoly0_E(:,2) = (/0.,6.,12.,17.,15./) call PPM_reconstruction(5, (/1.,1.,1.,1.,1./), (/0.,5.,7.,16.,15./), ppoly0_E(1:5,:), ppoly0_coefficients(1:5,:) ) - remapping_unit_tests = remapping_unit_tests .or. test_answer(5, ppoly0_E(:,1), (/0.,3.,6.,16.,15./), 'Limits PPM: left edges') - remapping_unit_tests = remapping_unit_tests .or. test_answer(5, ppoly0_E(:,2), (/0.,6.,9.,16.,15./), 'Limits PPM: right edges') - remapping_unit_tests = remapping_unit_tests .or. test_answer(5, ppoly0_coefficients(:,1), (/0.,3.,6.,16.,15./), 'Limits PPM: P0') - remapping_unit_tests = remapping_unit_tests .or. test_answer(5, ppoly0_coefficients(:,2), (/0.,6.,0.,0.,0./), 'Limits PPM: P1') - remapping_unit_tests = remapping_unit_tests .or. test_answer(5, ppoly0_coefficients(:,3), (/0.,-3.,3.,0.,0./), 'Limits PPM: P2') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 5, ppoly0_E(:,1), (/0.,3.,6.,16.,15./), 'Limits PPM: left edges') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 5, ppoly0_E(:,2), (/0.,6.,9.,16.,15./), 'Limits PPM: right edges') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 5, ppoly0_coefficients(:,1), (/0.,3.,6.,16.,15./), 'Limits PPM: P0') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 5, ppoly0_coefficients(:,2), (/0.,6.,0.,0.,0./), 'Limits PPM: P1') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 5, ppoly0_coefficients(:,3), (/0.,-3.,3.,0.,0./), 'Limits PPM: P2') + + call PLM_reconstruction(4, (/0.,1.,1.,0./), (/5.,4.,2.,1./), ppoly0_E(1:4,:), ppoly0_coefficients(1:4,:) ) + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 4, ppoly0_E(1:4,1), (/5.,5.,3.,1./), 'PPM: left edges h=0110') + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 4, ppoly0_E(1:4,2), (/5.,3.,1.,1./), 'PPM: right edges h=0110') + call remap_via_sub_cells( 4, (/0.,1.,1.,0./), (/5.,4.,2.,1./), ppoly0_E(1:4,:), ppoly0_coefficients(1:4,:), & + 2, (/1.,1./), INTEGRATION_PLM, .false., u2, err ) + remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 2, u2, (/4.,2./), 'PLM: remapped h=0110->h=11') deallocate(ppoly0_E, ppoly0_S, ppoly0_coefficients) - write(*,*) '==========================================================' + if (.not. remapping_unit_tests) write(*,*) 'Pass' - contains +end function remapping_unit_tests - !> Returns true if any cell of u and u_true are not identical. Returns false otherwise. - logical function test_answer(n, u, u_true, label) - integer, intent(in) :: n !< Number of cells in u - real, dimension(n), intent(in) :: u !< Values to test - real, dimension(n), intent(in) :: u_true !< Values to test against (correct answer) - character(len=*), intent(in) :: label !< Message - ! Local variables - integer :: k +!> Returns true if any cell of u and u_true are not identical. Returns false otherwise. +logical function test_answer(verbose, n, u, u_true, label) + logical, intent(in) :: verbose !< If true, write results to stdout + integer, intent(in) :: n !< Number of cells in u + real, dimension(n), intent(in) :: u !< Values to test + real, dimension(n), intent(in) :: u_true !< Values to test against (correct answer) + character(len=*), intent(in) :: label !< Message + ! Local variables + integer :: k - test_answer = .false. + test_answer = .false. + do k = 1, n + if (u(k) /= u_true(k)) test_answer = .true. + enddo + if (test_answer .or. verbose) then + write(*,'(a4,2a24,x,a)') 'k','Calculated value','Correct value',label do k = 1, n - if (u(k) /= u_true(k)) test_answer = .true. + if (u(k) /= u_true(k)) then + write(*,'(i4,1p2e24.16,a,1pe24.16,a)') k,u(k),u_true(k),' err=',u(k)-u_true(k),' < wrong' + else + write(*,'(i4,1p2e24.16)') k,u(k),u_true(k) + endif enddo - if (test_answer .or. .true.) then - write(*,'(a4,2a24,x,a)') 'k','Calculated value','Correct value',label - do k = 1, n - if (u(k) /= u_true(k)) then - write(*,'(i4,1p2e24.16,a,1pe24.16,a)') k,u(k),u_true(k),' err=',u(k)-u_true(k),' < wrong' - else - write(*,'(i4,1p2e24.16)') k,u(k),u_true(k) - endif - enddo - endif + endif - end function test_answer +end function test_answer - !> Convenience function for printing grid to screen - subroutine dumpGrid(n,h,x,u) +!> Convenience function for printing grid to screen +subroutine dumpGrid(n,h,x,u) integer, intent(in) :: n !< Number of cells real, dimension(:), intent(in) :: h !< Cell thickness real, dimension(:), intent(in) :: x !< Interface delta @@ -1765,8 +1778,6 @@ subroutine dumpGrid(n,h,x,u) write(*,'("i=",5x,20i10)') (i,i=1,n) write(*,'("h=",5x,20es10.2)') (h(i),i=1,n) write(*,'("u=",5x,20es10.2)') (u(i),i=1,n) - end subroutine dumpGrid - -end function remapping_unit_tests +end subroutine dumpGrid end module MOM_remapping diff --git a/src/ALE/coord_rho.F90 b/src/ALE/coord_rho.F90 index 15a08a77dc..b66f76b9cd 100644 --- a/src/ALE/coord_rho.F90 +++ b/src/ALE/coord_rho.F90 @@ -198,10 +198,10 @@ subroutine build_rho_column(CS, remapCS, nz, depth, h, T, S, eqn_of_state, zInte h1(k) = x1(k+1) - x1(k) end do - call remapping_core_h(nz, h0, S, nz, h1, Tmp, remapCS) + call remapping_core_h(remapCS, nz, h0, S, nz, h1, Tmp) S_tmp(:) = Tmp(:) - call remapping_core_h(nz, h0, T, nz, h1, Tmp, remapCS) + call remapping_core_h(remapCS, nz, h0, T, nz, h1, Tmp) T_tmp(:) = Tmp(:) ! Compute the deviation between two successive grids diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 2510aa1166..94e7fed7e8 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -60,6 +60,7 @@ module MOM use MOM_time_manager, only : time_type, set_time, time_type_to_real, operator(+) use MOM_time_manager, only : operator(-), operator(>), operator(*), operator(/) use MOM_time_manager, only : increment_date +use MOM_unit_tests, only : unit_tests ! MOM core modules use MOM_ALE, only : ALE_init, ALE_end, ALE_main, ALE_CS, adjustGridForIntegrity @@ -1790,6 +1791,13 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo "\t0 = Only FATAL messages\n" // & "\t2 = Only FATAL, WARNING, NOTE [default]\n" // & "\t9 = All)", default=2) + call get_param(param_file, "MOM", "DO_UNIT_TESTS", do_unit_tests, & + "If True, exercises unit tests at model start up.", & + default=.false.) + if (do_unit_tests) then + call unit_tests(verbosity) + endif + call get_param(param_file, "MOM", "SPLIT", CS%split, & "Use the split time stepping if true.", default=.true.) if (CS%split) then @@ -2535,11 +2543,6 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo .not.((dirs%input_filename(1:1) == 'r') .and. & (LEN_TRIM(dirs%input_filename) == 1)) - ! Undocumented parameter: set DO_UNIT_TESTS=True to invoke unit_tests s/r - ! which calls unit tests provided by some modules. - call get_param(param_file, "MOM", "DO_UNIT_TESTS", do_unit_tests, default=.false.) - if (do_unit_tests) call unit_tests - call callTree_leave("initialize_MOM()") call cpu_clock_end(id_clock_init) @@ -2585,32 +2588,6 @@ subroutine finish_MOM_initialization(Time, dirs, CS, fluxes) end subroutine finish_MOM_initialization - -!> Calls unit tests for other modules. These are NOT normally invoked -!! and so we provide the module use statements here rather than in the module -!! header. This is an exception to our usual coding standards. -!! Note that if a unit test returns true, a FATAL error is triggered. -subroutine unit_tests - use MOM_string_functions, only : string_functions_unit_tests - use MOM_remapping, only : remapping_unit_tests - use MOM_neutral_diffusion, only : neutral_diffusion_unit_tests - use MOM_diag_vkernels, only : diag_vkernels_unit_tests - - if (is_root_pe()) then ! The following need only be tested on 1 PE - if (string_functions_unit_tests()) call MOM_error(FATAL, & - "MOM/initialize_MOM/unit_tests: string_functions_unit_tests FAILED") - if (remapping_unit_tests()) call MOM_error(FATAL, & - "MOM/initialize_MOM/unit_tests: remapping_unit_tests FAILED") - if (neutral_diffusion_unit_tests()) call MOM_error(FATAL, & - "MOM/initialize_MOM/unit_tests: neutralDiffusionUnitTests FAILED") - if (diag_vkernels_unit_tests()) call MOM_error(FATAL, & - "MOM/initialize_MOM/unit_tests: diag_vkernels_unit_tests FAILED") - - endif - -end subroutine unit_tests - - !> Register the diagnostics subroutine register_diags(Time, G, GV, CS, ADp) type(time_type), intent(in) :: Time !< current model time diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index eedc0325d2..994725b8b9 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -2179,9 +2179,10 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time) ! Pretty sure we need to check for source/target grid consistency here segment%field(m)%buffer_dst(i,j,:)=0.0 ! initialize remap destination buffer if (G%mask2dT(i,j)>0.) then - call remapping_core_h(segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:),& - segment%field(m)%buffer_src(i,j,:),G%ke, h(i,j,:),& - segment%field(m)%buffer_dst(i,j,:),OBC%remap_CS) + call remapping_core_h(OBC%remap_CS, & + segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:), & + segment%field(m)%buffer_src(i,j,:), & + G%ke, h(i,j,:), segment%field(m)%buffer_dst(i,j,:)) endif enddo enddo diff --git a/src/core/MOM_unit_tests.F90 b/src/core/MOM_unit_tests.F90 new file mode 100644 index 0000000000..5417301c64 --- /dev/null +++ b/src/core/MOM_unit_tests.F90 @@ -0,0 +1,36 @@ +!> Invokes unit tests in all modules that have them +module MOM_unit_tests + +use MOM_error_handler, only : MOM_error, FATAL, is_root_pe + +use MOM_string_functions, only : string_functions_unit_tests +use MOM_remapping, only : remapping_unit_tests +use MOM_neutral_diffusion, only : neutral_diffusion_unit_tests +use MOM_diag_vkernels, only : diag_vkernels_unit_tests + +contains + +!> Calls unit tests for other modules. +!! Note that if a unit test returns true, a FATAL error is triggered. +subroutine unit_tests(verbosity) + ! Arguments + integer, intent(in) :: verbosity !< The verbosity level + ! Local variables + logical :: verbose + + verbose = verbosity>=5 + + if (is_root_pe()) then ! The following need only be tested on 1 PE + if (string_functions_unit_tests(verbose)) call MOM_error(FATAL, & + "MOM_unit_tests: string_functions_unit_tests FAILED") + if (remapping_unit_tests(verbose)) call MOM_error(FATAL, & + "MOM_unit_tests: remapping_unit_tests FAILED") + if (neutral_diffusion_unit_tests(verbose)) call MOM_error(FATAL, & + "MOM_unit_tests: neutralDiffusionUnitTests FAILED") + if (diag_vkernels_unit_tests(verbose)) call MOM_error(FATAL, & + "MOM_unit_tests: diag_vkernels_unit_tests FAILED") + endif + +end subroutine unit_tests + +end module MOM_unit_tests diff --git a/src/diagnostics/MOM_wave_speed.F90 b/src/diagnostics/MOM_wave_speed.F90 index 40c3c4281f..9b4f47e5ef 100644 --- a/src/diagnostics/MOM_wave_speed.F90 +++ b/src/diagnostics/MOM_wave_speed.F90 @@ -438,7 +438,7 @@ subroutine wave_speed(h, tv, G, GV, cg1, CS, full_halos, use_ebt_mode, & else mode_struct(1:kc)=0. endif - call remapping_core_h(kc, Hc, mode_struct, nz, h(i,j,:), modal_structure(i,j,:), CS%remapping_CS) + call remapping_core_h(CS%remapping_CS, kc, Hc, mode_struct, nz, h(i,j,:), modal_structure(i,j,:)) endif else cg1(i,j) = 0.0 diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index 66e5970d7f..d579150879 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -317,9 +317,8 @@ subroutine diag_remap_do_remap(remap_cs, G, h, staggered_in_x, staggered_in_y, & endif h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:)) h_dest(:) = 0.5 * (remap_cs%h(i,j,:) + remap_cs%h(i+1,j,:)) - call remapping_core_h(nz_src, h_src(:), field(I,j,:), & - nz_dest, h_dest(:), remapped_field(I,j,:), & - remap_cs%remap_cs) + call remapping_core_h(remap_cs%remap_cs, nz_src, h_src(:), field(I,j,:), & + nz_dest, h_dest(:), remapped_field(I,j,:)) if (mask_vanished_layers) then ! This only works for z-like output do k=1, nz_dest if (h_dest(k) == 0.) remapped_field(i, j, k:nz_dest) = missing_value @@ -336,9 +335,8 @@ subroutine diag_remap_do_remap(remap_cs, G, h, staggered_in_x, staggered_in_y, & endif h_src(:) = 0.5 * (h(i,j,:) + h(i,j+1,:)) h_dest(:) = 0.5 * (remap_cs%h(i,j,:) + remap_cs%h(i,j+1,:) ) - call remapping_core_h(nz_src, h_src(:), field(i,J,:), & - nz_dest, h_dest(:), remapped_field(i,J,:), & - remap_cs%remap_cs) + call remapping_core_h(remap_cs%remap_cs, nz_src, h_src(:), field(i,J,:), & + nz_dest, h_dest(:), remapped_field(i,J,:)) if (mask_vanished_layers) then ! This only works for z-like output do k=1, nz_dest if (h_dest(k) == 0.) remapped_field(i,j,k) = missing_value @@ -354,9 +352,8 @@ subroutine diag_remap_do_remap(remap_cs, G, h, staggered_in_x, staggered_in_y, & if (mask(i,j, 1) == 0.) cycle endif h_dest(:) = remap_cs%h(i,j,:) - call remapping_core_h(nz_src, h(i,j,:), field(i,j,:), & - nz_dest, h_dest(:), remapped_field(i,j,:), & - remap_cs%remap_cs) + call remapping_core_h(remap_cs%remap_cs, nz_src, h(i,j,:), field(i,j,:), & + nz_dest, h_dest(:), remapped_field(i,j,:)) if (mask_vanished_layers) then ! This only works for z-like output do k=1, nz_dest if (h_dest(k)==0.) remapped_field(i,j,k) = missing_value diff --git a/src/framework/MOM_diag_vkernels.F90 b/src/framework/MOM_diag_vkernels.F90 index 18d283bffc..1b8fb58b6d 100644 --- a/src/framework/MOM_diag_vkernels.F90 +++ b/src/framework/MOM_diag_vkernels.F90 @@ -10,9 +10,6 @@ module MOM_diag_vkernels public interpolate_column public reintegrate_column -! Module data used only for debugging -logical :: verbose_diag_vkernels_unit_tests = .true. ! If true, always write out unit tests results - contains !> Linearly interpolate interface data, u_src, from grid h_src to a grid h_dest @@ -168,187 +165,192 @@ subroutine reintegrate_column(nsrc, h_src, uh_src, ndest, h_dest, missing_value, end subroutine reintegrate_column !> Returns true if any unit tests for module MOM_diag_vkernels fail -logical function diag_vkernels_unit_tests() +logical function diag_vkernels_unit_tests(verbose) + logical, intent(in) :: verbose !< If true, write results to stdout ! Local variables - real, parameter :: missing_value=-9.999999999E9 ! Value to use for vanished layers - logical :: fail + real, parameter :: mv=-9.999999999E9 ! Value to use for vanished layers + logical :: fail, v + + v = verbose - write(0,*) '============ MOM_diag_kernels: diag_vkernels_unit_tests ==================' - write(0,*) '- - - - - - - - - - interpolation tests - - - - - - - - - - - - - - - - -' + write(0,*) '==== MOM_diag_kernels: diag_vkernels_unit_tests ==========' + if (v) write(0,*) '- - - - - - - - - - interpolation tests - - - - - - - - -' - fail = test_interp('Identity: 3 layer', & + fail = test_interp(v,mv,'Identity: 3 layer', & 3, (/1.,2.,3./), (/1.,2.,3.,4./), & 3, (/1.,2.,3./), (/1.,2.,3.,4./) ) diag_vkernels_unit_tests = fail - fail = test_interp('A: 3 layer to 2', & + fail = test_interp(v,mv,'A: 3 layer to 2', & 3, (/1.,1.,1./), (/1.,2.,3.,4./), & 2, (/1.5,1.5/), (/1.,2.5,4./) ) diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail - fail = test_interp('B: 2 layer to 3', & + fail = test_interp(v,mv,'B: 2 layer to 3', & 2, (/1.5,1.5/), (/1.,4.,7./), & 3, (/1.,1.,1./), (/1.,3.,5.,7./) ) diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail - fail = test_interp('C: 3 layer (vanished middle) to 2', & + fail = test_interp(v,mv,'C: 3 layer (vanished middle) to 2', & 3, (/1.,0.,2./), (/1.,2.,2.,3./), & 2, (/1.,2./), (/1.,2.,3./) ) diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail - fail = test_interp('D: 3 layer (deep) to 3', & + fail = test_interp(v,mv,'D: 3 layer (deep) to 3', & 3, (/1.,2.,3./), (/1.,2.,4.,7./), & 2, (/2.,2./), (/1.,3.,5./) ) diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail - fail = test_interp('E: 3 layer to 3 (deep)', & + fail = test_interp(v,mv,'E: 3 layer to 3 (deep)', & 3, (/1.,2.,4./), (/1.,2.,4.,8./), & 3, (/2.,3.,4./), (/1.,3.,6.,8./) ) diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail - fail = test_interp('F: 3 layer to 4 with vanished top/botton', & + fail = test_interp(v,mv,'F: 3 layer to 4 with vanished top/botton', & 3, (/1.,2.,4./), (/1.,2.,4.,8./), & - 4, (/0.,2.,5.,0./), (/missing_value,1.,3.,8.,missing_value/) ) + 4, (/0.,2.,5.,0./), (/mv,1.,3.,8.,mv/) ) diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail - fail = test_interp('Fs: 3 layer to 4 with vanished top/botton (shallow)', & + fail = test_interp(v,mv,'Fs: 3 layer to 4 with vanished top/botton (shallow)', & 3, (/1.,2.,4./), (/1.,2.,4.,8./), & - 4, (/0.,2.,4.,0./), (/missing_value,1.,3.,7.,missing_value/) ) + 4, (/0.,2.,4.,0./), (/mv,1.,3.,7.,mv/) ) diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail - fail = test_interp('Fd: 3 layer to 4 with vanished top/botton (deep)', & + fail = test_interp(v,mv,'Fd: 3 layer to 4 with vanished top/botton (deep)', & 3, (/1.,2.,4./), (/1.,2.,4.,8./), & - 4, (/0.,2.,6.,0./), (/missing_value,1.,3.,8.,missing_value/) ) + 4, (/0.,2.,6.,0./), (/mv,1.,3.,8.,mv/) ) diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail - write(0,*) '- - - - - - - - - - reintegration tests - - - - - - - - - - - - - - - - -' + if (v) write(0,*) '- - - - - - - - - - reintegration tests - - - - - - - - -' - fail = test_reintegrate('Identity: 3 layer', & + fail = test_reintegrate(v,mv,'Identity: 3 layer', & 3, (/1.,2.,3./), (/-5.,2.,1./), & 3, (/1.,2.,3./), (/-5.,2.,1./) ) diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail - fail = test_reintegrate('A: 3 layer to 2', & + fail = test_reintegrate(v,mv,'A: 3 layer to 2', & 3, (/2.,2.,2./), (/-5.,2.,1./), & 2, (/3.,3./), (/-4.,2./) ) diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail - fail = test_reintegrate('A: 3 layer to 2 (deep)', & + fail = test_reintegrate(v,mv,'A: 3 layer to 2 (deep)', & 3, (/2.,2.,2./), (/-5.,2.,1./), & 2, (/3.,4./), (/-4.,2./) ) diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail - fail = test_reintegrate('A: 3 layer to 2 (shallow)', & + fail = test_reintegrate(v,mv,'A: 3 layer to 2 (shallow)', & 3, (/2.,2.,2./), (/-5.,2.,1./), & 2, (/3.,2./), (/-4.,1.5/) ) diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail - fail = test_reintegrate('B: 3 layer to 4 with vanished top/bottom', & + fail = test_reintegrate(v,mv,'B: 3 layer to 4 with vanished top/bottom', & 3, (/2.,2.,2./), (/-5.,2.,1./), & 4, (/0.,3.,3.,0./), (/0.,-4.,2.,0./) ) diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail - fail = test_reintegrate('C: 3 layer to 4 with vanished top//middle/bottom', & + fail = test_reintegrate(v,mv,'C: 3 layer to 4 with vanished top//middle/bottom', & 3, (/2.,2.,2./), (/-5.,2.,1./), & 5, (/0.,3.,0.,3.,0./), (/0.,-4.,0.,2.,0./) ) diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail - fail = test_reintegrate('D: 3 layer to 3 (vanished)', & + fail = test_reintegrate(v,mv,'D: 3 layer to 3 (vanished)', & 3, (/2.,2.,2./), (/-5.,2.,1./), & 3, (/0.,0.,0./), (/0.,0.,0./) ) diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail - fail = test_reintegrate('D: 3 layer (vanished) to 3', & + fail = test_reintegrate(v,mv,'D: 3 layer (vanished) to 3', & 3, (/0.,0.,0./), (/-5.,2.,1./), & - 3, (/2.,2.,2./), (/missing_value, missing_value, missing_value/) ) + 3, (/2.,2.,2./), (/mv, mv, mv/) ) diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail - fail = test_reintegrate('D: 3 layer (vanished) to 3 (vanished)', & + fail = test_reintegrate(v,mv,'D: 3 layer (vanished) to 3 (vanished)', & 3, (/0.,0.,0./), (/-5.,2.,1./), & - 3, (/0.,0.,0./), (/missing_value, missing_value, missing_value/) ) + 3, (/0.,0.,0./), (/mv, mv, mv/) ) diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail - fail = test_reintegrate('D: 3 layer (vanished) to 3 (vanished)', & + fail = test_reintegrate(v,mv,'D: 3 layer (vanished) to 3 (vanished)', & 3, (/0.,0.,0./), (/0.,0.,0./), & - 3, (/0.,0.,0./), (/missing_value, missing_value, missing_value/) ) + 3, (/0.,0.,0./), (/mv, mv, mv/) ) diag_vkernels_unit_tests = diag_vkernels_unit_tests .or. fail - write(0,*) '==========================================================================' - - contains - - !> Returns true if a test of interpolate_column() produces the wrong answer - logical function test_interp(msg, nsrc, h_src, u_src, ndest, h_dest, u_true) - character(len=*), intent(in) :: msg !< Message to label test - integer, intent(in) :: nsrc !< Number of source cells - real, dimension(nsrc), intent(in) :: h_src !< Thickness of source cells - real, dimension(nsrc+1), intent(in) :: u_src !< Values at source cell interfaces - integer, intent(in) :: ndest !< Number of destination cells - real, dimension(ndest), intent(in) :: h_dest !< Thickness of destination cells - real, dimension(ndest+1), intent(in) :: u_true !< Correct value at destination cell interfaces - ! Local variables - real, dimension(ndest+1) :: u_dest ! Interpolated value at destination cell interfaces - integer :: k - real :: error - logical :: print_results - - ! Interpolate from src to dest - call interpolate_column(nsrc, h_src, u_src, ndest, h_dest, missing_value, u_dest) - - test_interp = .false. + if (.not. fail) write(*,*) 'Pass' + +end function diag_vkernels_unit_tests + +!> Returns true if a test of interpolate_column() produces the wrong answer +logical function test_interp(verbose, missing_value, msg, nsrc, h_src, u_src, ndest, h_dest, u_true) + logical, intent(in) :: verbose !< If true, write results to stdout + real, intent(in) :: missing_value !< Value to indicate missing data + character(len=*), intent(in) :: msg !< Message to label test + integer, intent(in) :: nsrc !< Number of source cells + real, dimension(nsrc), intent(in) :: h_src !< Thickness of source cells + real, dimension(nsrc+1), intent(in) :: u_src !< Values at source cell interfaces + integer, intent(in) :: ndest !< Number of destination cells + real, dimension(ndest), intent(in) :: h_dest !< Thickness of destination cells + real, dimension(ndest+1), intent(in) :: u_true !< Correct value at destination cell interfaces + ! Local variables + real, dimension(ndest+1) :: u_dest ! Interpolated value at destination cell interfaces + integer :: k + real :: error + logical :: print_results + + ! Interpolate from src to dest + call interpolate_column(nsrc, h_src, u_src, ndest, h_dest, missing_value, u_dest) + + test_interp = .false. + do k=1,ndest+1 + if (u_dest(k)/=u_true(k)) test_interp = .true. + enddo + if (verbose .or. test_interp) then + write(0,'(2a)') ' Test: ',msg + write(0,'(a3,3(a24))') 'k','u_result','u_true','error' do k=1,ndest+1 - if (u_dest(k)/=u_true(k)) test_interp = .true. + error = u_dest(k)-u_true(k) + if (error==0.) then + write(0,'(i3,3(1pe24.16))') k,u_dest(k),u_true(k),u_dest(k)-u_true(k) + else + write(0,'(i3,3(1pe24.16),x,a)') k,u_dest(k),u_true(k),u_dest(k)-u_true(k),'<--- WRONG!' + endif enddo - if (verbose_diag_vkernels_unit_tests .or. test_interp) then - write(0,'(2a)') ' Test: ',msg - write(0,'(a3,3(a24))') 'k','u_result','u_true','error' - do k=1,ndest+1 - error = u_dest(k)-u_true(k) - if (error==0.) then - write(0,'(i3,3(1pe24.16))') k,u_dest(k),u_true(k),u_dest(k)-u_true(k) - else - write(0,'(i3,3(1pe24.16),x,a)') k,u_dest(k),u_true(k),u_dest(k)-u_true(k),'<--- WRONG!' - endif - enddo - endif - end function test_interp - - !> Returns true if a test of reintegrate_column() produces the wrong answer - logical function test_reintegrate(msg, nsrc, h_src, uh_src, ndest, h_dest, uh_true) - character(len=*), intent(in) :: msg !< Message to label test - integer, intent(in) :: nsrc !< Number of source cells - real, dimension(nsrc), intent(in) :: h_src !< Thickness of source cells - real, dimension(nsrc), intent(in) :: uh_src !< Values of source cell stuff - integer, intent(in) :: ndest !< Number of destination cells - real, dimension(ndest), intent(in) :: h_dest !< Thickness of destination cells - real, dimension(ndest), intent(in) :: uh_true !< Correct value of destination cell stuff - ! Local variables - real, dimension(ndest) :: uh_dest ! Reintegrated value on destination cells - integer :: k - real :: error - logical :: print_results - - ! Interpolate from src to dest - call reintegrate_column(nsrc, h_src, uh_src, ndest, h_dest, missing_value, uh_dest) - - test_reintegrate = .false. + endif +end function test_interp + +!> Returns true if a test of reintegrate_column() produces the wrong answer +logical function test_reintegrate(verbose, missing_value, msg, nsrc, h_src, uh_src, ndest, h_dest, uh_true) + logical, intent(in) :: verbose !< If true, write results to stdout + real, intent(in) :: missing_value !< Value to indicate missing data + character(len=*), intent(in) :: msg !< Message to label test + integer, intent(in) :: nsrc !< Number of source cells + real, dimension(nsrc), intent(in) :: h_src !< Thickness of source cells + real, dimension(nsrc), intent(in) :: uh_src !< Values of source cell stuff + integer, intent(in) :: ndest !< Number of destination cells + real, dimension(ndest), intent(in) :: h_dest !< Thickness of destination cells + real, dimension(ndest), intent(in) :: uh_true !< Correct value of destination cell stuff + ! Local variables + real, dimension(ndest) :: uh_dest ! Reintegrated value on destination cells + integer :: k + real :: error + logical :: print_results + + ! Interpolate from src to dest + call reintegrate_column(nsrc, h_src, uh_src, ndest, h_dest, missing_value, uh_dest) + + test_reintegrate = .false. + do k=1,ndest + if (uh_dest(k)/=uh_true(k)) test_reintegrate = .true. + enddo + if (verbose .or. test_reintegrate) then + write(0,'(2a)') ' Test: ',msg + write(0,'(a3,3(a24))') 'k','uh_result','uh_true','error' do k=1,ndest - if (uh_dest(k)/=uh_true(k)) test_reintegrate = .true. + error = uh_dest(k)-uh_true(k) + if (error==0.) then + write(0,'(i3,3(1pe24.16))') k,uh_dest(k),uh_true(k),uh_dest(k)-uh_true(k) + else + write(0,'(i3,3(1pe24.16),x,a)') k,uh_dest(k),uh_true(k),uh_dest(k)-uh_true(k),'<--- WRONG!' + endif enddo - if (verbose_diag_vkernels_unit_tests .or. test_reintegrate) then - write(0,'(2a)') ' Test: ',msg - write(0,'(a3,3(a24))') 'k','uh_result','uh_true','error' - do k=1,ndest - error = uh_dest(k)-uh_true(k) - if (error==0.) then - write(0,'(i3,3(1pe24.16))') k,uh_dest(k),uh_true(k),uh_dest(k)-uh_true(k) - else - write(0,'(i3,3(1pe24.16),x,a)') k,uh_dest(k),uh_true(k),uh_dest(k)-uh_true(k),'<--- WRONG!' - endif - enddo - endif - end function test_reintegrate - -end function diag_vkernels_unit_tests + endif +end function test_reintegrate end module MOM_diag_vkernels diff --git a/src/framework/MOM_string_functions.F90 b/src/framework/MOM_string_functions.F90 index 63b5063d5a..46431df09b 100644 --- a/src/framework/MOM_string_functions.F90 +++ b/src/framework/MOM_string_functions.F90 @@ -331,75 +331,88 @@ end function extract_real end function remove_spaces !> Returns true if a unit test of string_functions fails. -logical function string_functions_unit_tests() +logical function string_functions_unit_tests(verbose) + ! Arguments + logical, intent(in) :: verbose !< If true, write results to stdout + ! Local variables integer :: i(5) = (/ -1, 1, 3, 3, 0 /) real :: r(8) = (/ 0., 1., -2., 1.3, 3.E-11, 3.E-11, 3.E-11, -5.1E12 /) - logical :: fail, verbose - verbose = .false. + logical :: fail, v fail = .false. + v = verbose write(*,*) '==== MOM_string_functions: string_functions_unit_tests ===' - fail = fail .or. localTestS(left_int(-1),'-1') - fail = fail .or. localTestS(left_ints(i(:)),'-1, 1, 3, 3, 0') - fail = fail .or. localTestS(left_real(0.),'0.0') - fail = fail .or. localTestS(left_reals(r(:)),'0.0, 1.0, -2.0, 1.3, 3*3.0E-11, -5.1E+12') - fail = fail .or. localTestS(left_reals(r(:),sep=' '),'0.0 1.0 -2.0 1.3 3*3.0E-11 -5.1E+12') - fail = fail .or. localTestS(left_reals(r(:),sep=','),'0.0,1.0,-2.0,1.3,3*3.0E-11,-5.1E+12') - fail = fail .or. localTestS(extractWord("One Two,Three",1),"One") - fail = fail .or. localTestS(extractWord("One Two,Three",2),"Two") - fail = fail .or. localTestS(extractWord("One Two,Three",3),"Three") - fail = fail .or. localTestS(extractWord("One Two, Three",3),"Three") - fail = fail .or. localTestS(extractWord(" One Two,Three",1),"One") - fail = fail .or. localTestS(extract_word("One,Two,Three",",",3),"Three") - fail = fail .or. localTestS(extract_word("One,Two,Three",",",4),"") - fail = fail .or. localTestS(remove_spaces("1 2 3"),"123") - fail = fail .or. localTestS(remove_spaces(" 1 2 3"),"123") - fail = fail .or. localTestS(remove_spaces("1 2 3 "),"123") - fail = fail .or. localTestS(remove_spaces("123"),"123") - fail = fail .or. localTestS(remove_spaces(" "),"") - fail = fail .or. localTestS(remove_spaces(""),"") - fail = fail .or. localTestI(extract_integer("1","",1),1) - fail = fail .or. localTestI(extract_integer("1,2,3",",",1),1) - fail = fail .or. localTestI(extract_integer("1,2",",",2),2) - fail = fail .or. localTestI(extract_integer("1,2",",",3),0) - fail = fail .or. localTestI(extract_integer("1,2",",",4,4),4) - fail = fail .or. localTestR(extract_real("1.","",1),1.) - fail = fail .or. localTestR(extract_real("1.,2.,3.",",",1),1.) - fail = fail .or. localTestR(extract_real("1.,2.",",",2),2.) - fail = fail .or. localTestR(extract_real("1.,2.",",",3),0.) - fail = fail .or. localTestR(extract_real("1.,2.",",",4,4.),4.) + fail = fail .or. localTestS(v,left_int(-1),'-1') + fail = fail .or. localTestS(v,left_ints(i(:)),'-1, 1, 3, 3, 0') + fail = fail .or. localTestS(v,left_real(0.),'0.0') + fail = fail .or. localTestS(v,left_reals(r(:)),'0.0, 1.0, -2.0, 1.3, 3*3.0E-11, -5.1E+12') + fail = fail .or. localTestS(v,left_reals(r(:),sep=' '),'0.0 1.0 -2.0 1.3 3*3.0E-11 -5.1E+12') + fail = fail .or. localTestS(v,left_reals(r(:),sep=','),'0.0,1.0,-2.0,1.3,3*3.0E-11,-5.1E+12') + fail = fail .or. localTestS(v,extractWord("One Two,Three",1),"One") + fail = fail .or. localTestS(v,extractWord("One Two,Three",2),"Two") + fail = fail .or. localTestS(v,extractWord("One Two,Three",3),"Three") + fail = fail .or. localTestS(v,extractWord("One Two, Three",3),"Three") + fail = fail .or. localTestS(v,extractWord(" One Two,Three",1),"One") + fail = fail .or. localTestS(v,extract_word("One,Two,Three",",",3),"Three") + fail = fail .or. localTestS(v,extract_word("One,Two,Three",",",4),"") + fail = fail .or. localTestS(v,remove_spaces("1 2 3"),"123") + fail = fail .or. localTestS(v,remove_spaces(" 1 2 3"),"123") + fail = fail .or. localTestS(v,remove_spaces("1 2 3 "),"123") + fail = fail .or. localTestS(v,remove_spaces("123"),"123") + fail = fail .or. localTestS(v,remove_spaces(" "),"") + fail = fail .or. localTestS(v,remove_spaces(""),"") + fail = fail .or. localTestI(v,extract_integer("1","",1),1) + fail = fail .or. localTestI(v,extract_integer("1,2,3",",",1),1) + fail = fail .or. localTestI(v,extract_integer("1,2",",",2),2) + fail = fail .or. localTestI(v,extract_integer("1,2",",",3),0) + fail = fail .or. localTestI(v,extract_integer("1,2",",",4,4),4) + fail = fail .or. localTestR(v,extract_real("1.","",1),1.) + fail = fail .or. localTestR(v,extract_real("1.,2.,3.",",",1),1.) + fail = fail .or. localTestR(v,extract_real("1.,2.",",",2),2.) + fail = fail .or. localTestR(v,extract_real("1.,2.",",",3),0.) + fail = fail .or. localTestR(v,extract_real("1.,2.",",",4,4.),4.) if (.not. fail) write(*,*) 'Pass' - write(*,*) '==========================================================' string_functions_unit_tests = fail - contains - logical function localTestS(str1,str2) - character(len=*) :: str1, str2 - localTestS=.false. - if (trim(str1)/=trim(str2)) localTestS=.true. - if (localTestS .or. verbose) then - write(*,*) '>'//trim(str1)//'<' - if (localTestS) write(*,*) trim(str1),':',trim(str2), '<-- FAIL' - endif - end function localTestS - logical function localTestI(i1,i2) - integer :: i1,i2 - localTestI=.false. - if (i1/=i2) localTestI=.true. - if (localTestI .or. verbose) then - write(*,*) i1,i2 - if (localTestI) write(*,*) i1,'!=',i2, '<-- FAIL' - endif - end function localTestI - logical function localTestR(r1,r2) - real :: r1,r2 - localTestR=.false. - if (r1/=r2) localTestR=.true. - if (localTestR .or. verbose) then - write(*,*) r1,r2 - if (localTestR) write(*,*) r1,'!=',r2, '<-- FAIL' - endif - end function localTestR end function string_functions_unit_tests +!> True if str1 does not match str2. False otherwise. +logical function localTestS(verbose,str1,str2) + logical, intent(in) :: verbose !< If true, write results to stdout + character(len=*), intent(in) :: str1 !< String + character(len=*), intent(in) :: str2 !< String + localTestS=.false. + if (trim(str1)/=trim(str2)) localTestS=.true. + if (localTestS .or. verbose) then + write(*,*) '>'//trim(str1)//'<' + if (localTestS) write(*,*) trim(str1),':',trim(str2), '<-- FAIL' + endif +end function localTestS + +!> True if i1 is not equal to i2. False otherwise. +logical function localTestI(verbose,i1,i2) + logical, intent(in) :: verbose !< If true, write results to stdout + integer, intent(in) :: i1 !< Integer + integer, intent(in) :: i2 !< Integer + localTestI=.false. + if (i1/=i2) localTestI=.true. + if (localTestI .or. verbose) then + write(*,*) i1,i2 + if (localTestI) write(*,*) i1,'!=',i2, '<-- FAIL' + endif +end function localTestI + +!> True if r1 is not equal to r2. False otherwise. +logical function localTestR(verbose,r1,r2) + logical, intent(in) :: verbose !< If true, write results to stdout + real, intent(in) :: r1 !< Float + real, intent(in) :: r2 !< Float + localTestR=.false. + if (r1/=r2) localTestR=.true. + if (localTestR .or. verbose) then + write(*,*) r1,r2 + if (localTestR) write(*,*) r1,'!=',r2, '<-- FAIL' + endif +end function localTestR + !> Returns a directory name that is terminated with a "/" or "./" if the !! argument is an empty string. function slasher(dir) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index cc33ddae1d..24b409c5bb 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -997,8 +997,8 @@ subroutine cut_off_column_top(nk, tv, Rho0, G_earth, depth, min_thickness, & T0(k) = T(nk+1-k) h1(k) = h(nk+1-k) enddo - call remapping_core_h(nk, h0, T0, nk, h1, T1, remap_CS ) - call remapping_core_h(nk, h0, S0, nk, h1, S1, remap_CS ) + call remapping_core_h(remap_CS, nk, h0, T0, nk, h1, T1) + call remapping_core_h(remap_CS, nk, h0, S0, nk, h1, S1) do k=1,nk S(k) = S1(nk+1-k) T(k) = T1(nk+1-k) diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index 096a0cc58f..dc66b11861 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -382,8 +382,9 @@ subroutine apply_ALE_sponge(h, dt, G, CS) I1pdamp = 1.0 / (1.0 + damp) do m=1,CS%fldno tmp_val2(1:CS%nz_data) = CS%Ref_val(m)%p(1:CS%nz_data,c) - call remapping_core_h(CS%nz_data, CS%Ref_h(:,c), tmp_val2, & - CS%nz, h(i,j,:), tmp_val1, CS%remap_cs) + call remapping_core_h(CS%remap_cs, & + CS%nz_data, CS%Ref_h(:,c), tmp_val2, & + CS%nz, h(i,j,:), tmp_val1) !Backward Euler method CS%var(m)%p(i,j,:) = I1pdamp * & @@ -412,8 +413,9 @@ subroutine apply_ALE_sponge(h, dt, G, CS) damp = dt*CS%Iresttime_col_u(c) I1pdamp = 1.0 / (1.0 + damp) tmp_val2(1:CS%nz_data) = CS%Ref_val_u(1:CS%nz_data,c) - call remapping_core_h(CS%nz_data, CS%Ref_hu(:,c), tmp_val2, & - CS%nz, hu(i,j,:), tmp_val1, CS%remap_cs) + call remapping_core_h(CS%remap_cs, & + CS%nz_data, CS%Ref_hu(:,c), tmp_val2, & + CS%nz, hu(i,j,:), tmp_val1) !Backward Euler method CS%var_u(i,j,:) = I1pdamp * (CS%var_u(i,j,:) + tmp_val1 * damp) @@ -429,8 +431,9 @@ subroutine apply_ALE_sponge(h, dt, G, CS) damp = dt*CS%Iresttime_col_v(c) I1pdamp = 1.0 / (1.0 + damp) tmp_val2(1:CS%nz_data) = CS%Ref_val_v(1:CS%nz_data,c) - call remapping_core_h(CS%nz_data, CS%Ref_hv(:,c), tmp_val2, & - CS%nz, hv(i,j,:), tmp_val1, CS%remap_cs) + call remapping_core_h(CS%remap_cs, & + CS%nz_data, CS%Ref_hv(:,c), tmp_val2, & + CS%nz, hv(i,j,:), tmp_val1) !Backward Euler method CS%var_v(i,j,:) = I1pdamp * (CS%var_v(i,j,:) + tmp_val1 * damp) enddo diff --git a/src/tracer/MOM_neutral_diffusion.F90 b/src/tracer/MOM_neutral_diffusion.F90 index c039d5942f..7fb4abd307 100644 --- a/src/tracer/MOM_neutral_diffusion.F90 +++ b/src/tracer/MOM_neutral_diffusion.F90 @@ -9,7 +9,6 @@ module MOM_neutral_diffusion use MOM_diag_mediator, only : post_data, register_diag_field use MOM_EOS, only : EOS_type, calculate_compress, calculate_density_derivs use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe -use MOM_error_handler, only : MOM_get_verbosity use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_file_parser, only : openParameterBlock, closeParameterBlock use MOM_grid, only : ocean_grid_type @@ -1046,7 +1045,9 @@ subroutine neutral_surface_flux(nk, hl, hr, Tl, Tr, PiL, PiR, KoL, KoR, hEff, Fl end subroutine neutral_surface_flux !> Returns true if unit tests of neutral_diffusion functions fail. Otherwise returns false. -logical function neutral_diffusion_unit_tests() +logical function neutral_diffusion_unit_tests(verbose) + logical, intent(in) :: verbose !< It true, write results to stdout + ! Local variables integer, parameter :: nk = 4 real, dimension(nk+1) :: TiL, TiR1, TiR2, TiR4, Tio ! Test interface temperatures real, dimension(nk) :: TL ! Test layer temperatures @@ -1056,49 +1057,49 @@ logical function neutral_diffusion_unit_tests() integer, dimension(2*nk+2) :: KoL, KoR ! Test indexes real, dimension(2*nk+1) :: hEff ! Test positions real, dimension(2*nk+1) :: Flx ! Test flux + integer :: k + logical :: v - integer :: k, verbosity - - verbosity = MOM_get_verbosity() + v = verbose neutral_diffusion_unit_tests = .false. ! Normally return false - write(*,'(a)') '===== MOM_neutral_diffusion: neutral_diffusion_unit_tests ===============' - - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(1.,1.,1., 0.,1.,2., 1., 'FV: Straight line on uniform grid') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(1.,1.,0., 0.,4.,8., 7., 'FV: Vanished right cell') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(0.,1.,1., 0.,4.,8., 7., 'FV: Vanished left cell') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(1.,2.,4., 0.,3.,9., 4., 'FV: Stretched grid') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(2.,0.,2., 0.,1.,2., 0., 'FV: Vanished middle cell') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(0.,1.,0., 0.,1.,2., 2., 'FV: Vanished on both sides') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(1.,0.,0., 0.,1.,2., 0., 'FV: Two vanished cell sides') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(0.,0.,0., 0.,1.,2., 0., 'FV: All vanished cells') - - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(1.,1.,1., 0.,1.,2., 1., 'LSQ: Straight line on uniform grid') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(1.,1.,0., 0.,1.,2., 1., 'LSQ: Vanished right cell') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(0.,1.,1., 0.,1.,2., 1., 'LSQ: Vanished left cell') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(1.,2.,4., 0.,3.,9., 2., 'LSQ: Stretched grid') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(1.,0.,1., 0.,1.,2., 2., 'LSQ: Vanished middle cell') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(0.,1.,0., 0.,1.,2., 0., 'LSQ: Vanished on both sides') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(1.,0.,0., 0.,1.,2., 0., 'LSQ: Two vanished cell sides') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(0.,0.,0., 0.,1.,2., 0., 'LSQ: All vanished cells') + write(*,*) '==== MOM_neutral_diffusion: neutral_diffusion_unit_tests =' + + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(v,1.,1.,1., 0.,1.,2., 1., 'FV: Straight line on uniform grid') + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(v,1.,1.,0., 0.,4.,8., 7., 'FV: Vanished right cell') + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(v,0.,1.,1., 0.,4.,8., 7., 'FV: Vanished left cell') + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(v,1.,2.,4., 0.,3.,9., 4., 'FV: Stretched grid') + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(v,2.,0.,2., 0.,1.,2., 0., 'FV: Vanished middle cell') + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(v,0.,1.,0., 0.,1.,2., 2., 'FV: Vanished on both sides') + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(v,1.,0.,0., 0.,1.,2., 0., 'FV: Two vanished cell sides') + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(v,0.,0.,0., 0.,1.,2., 0., 'FV: All vanished cells') + + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(v,1.,1.,1., 0.,1.,2., 1., 'LSQ: Straight line on uniform grid') + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(v,1.,1.,0., 0.,1.,2., 1., 'LSQ: Vanished right cell') + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(v,0.,1.,1., 0.,1.,2., 1., 'LSQ: Vanished left cell') + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(v,1.,2.,4., 0.,3.,9., 2., 'LSQ: Stretched grid') + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(v,1.,0.,1., 0.,1.,2., 2., 'LSQ: Vanished middle cell') + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(v,0.,1.,0., 0.,1.,2., 0., 'LSQ: Vanished on both sides') + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(v,1.,0.,0., 0.,1.,2., 0., 'LSQ: Two vanished cell sides') + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(v,0.,0.,0., 0.,1.,2., 0., 'LSQ: All vanished cells') call interface_scalar(4, (/10.,10.,10.,10./), (/24.,18.,12.,6./), Tio, 1) !neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(5, Tio, (/27.,21.,15.,9.,3./), 'Linear profile, interface temperatures') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(5, Tio, (/24.,22.5,15.,7.5,6./), 'Linear profile, linear interface temperatures') + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(v,5, Tio, (/24.,22.5,15.,7.5,6./), 'Linear profile, linear interface temperatures') call interface_scalar(4, (/10.,10.,10.,10./), (/24.,18.,12.,6./), Tio, 2) - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(5, Tio, (/24.,22.,15.,8.,6./), 'Linear profile, PPM interface temperatures') - - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(-1.0, 0., 1.0, 1.0, 0.5, 'Check mid-point') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp( 0.0, 0., 1.0, 1.0, 0.0, 'Check bottom') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp( 0.1, 0., 1.1, 1.0, 0.0, 'Check below') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(-1.0, 0., 0.0, 1.0, 1.0, 'Check top') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(-1.0, 0., -0.1, 1.0, 1.0, 'Check above') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(-1.0, 0., 3.0, 1.0, 0.25, 'Check 1/4') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(-3.0, 0., 1.0, 1.0, 0.75, 'Check 3/4') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp( 1.0, 0., 1.0, 1.0, 0.0, 'Check dRho=0 below') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(-1.0, 0., -1.0, 1.0, 1.0, 'Check dRho=0 above') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp( 0.0, 0., 0.0, 1.0, 0.5, 'Check dRho=0 mid') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(-2.0, .5, 5.0, 0.5, 0.5, 'Check dP=0') + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(v,5, Tio, (/24.,22.,15.,8.,6./), 'Linear profile, PPM interface temperatures') + + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v,-1.0, 0., 1.0, 1.0, 0.5, 'Check mid-point') + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v, 0.0, 0., 1.0, 1.0, 0.0, 'Check bottom') + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v, 0.1, 0., 1.1, 1.0, 0.0, 'Check below') + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v,-1.0, 0., 0.0, 1.0, 1.0, 'Check top') + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v,-1.0, 0., -0.1, 1.0, 1.0, 'Check above') + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v,-1.0, 0., 3.0, 1.0, 0.25, 'Check 1/4') + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v,-3.0, 0., 1.0, 1.0, 0.75, 'Check 3/4') + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v, 1.0, 0., 1.0, 1.0, 0.0, 'Check dRho=0 below') + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v,-1.0, 0., -1.0, 1.0, 1.0, 'Check dRho=0 above') + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v, 0.0, 0., 0.0, 1.0, 0.5, 'Check dRho=0 mid') + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v,-2.0, .5, 5.0, 0.5, 0.5, 'Check dP=0') ! Identical columns call find_neutral_surface_positions(3, & @@ -1107,28 +1108,28 @@ logical function neutral_diffusion_unit_tests() (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Right positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS PiLRo, PiRLo, KoL, KoR, hEff) - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(3, KoL, KoR, PiLRo, PiRLo, hEff, & + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, KoL, KoR, PiLRo, PiRLo, hEff, & (/1,1,2,2,3,3,3,3/), & ! KoL (/1,1,2,2,3,3,3,3/), & ! KoR (/0.,0.,0.,0.,0.,0.,1.,1./), & ! pL (/0.,0.,0.,0.,0.,0.,1.,1./), & ! pR (/0.,10.,0.,10.,0.,10.,0./), & ! hEff 'Identical columns') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(8, & + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(v, 8, & absolute_positions(3, (/0.,10.,20.,30./), KoL, PiLRo), & (/0.,0.,10.,10.,20.,20.,30.,30./), '... left positions') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(8, & + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(v, 8, & absolute_positions(3, (/0.,10.,20.,30./), KoR, PiRLo), & (/0.,0.,10.,10.,20.,20.,30.,30./), '... right positions') call neutral_surface_flux(3, (/10.,10.,10./), (/10.,10.,10./), & ! nk, hL, hR (/20.,16.,12./), (/20.,16.,12./), & ! Tl, Tr PiLRo, PiRLo, KoL, KoR, hEff, Flx) - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(7, Flx, & + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(v, 7, Flx, & (/0.,0.,0.,0.,0.,0.,0./), 'Identical columns, rho flux (=0)') call neutral_surface_flux(3, (/10.,10.,10./), (/10.,10.,10./), & ! nk, hL, hR (/-1.,-1.,-1./), (/1.,1.,1./), & ! Sl, Sr PiLRo, PiRLo, KoL, KoR, hEff, Flx) - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(7, Flx, & + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(v, 7, Flx, & (/0.,20.,0.,20.,0.,20.,0./), 'Identical columns, S flux') ! Right column slightly cooler than left @@ -1138,17 +1139,17 @@ logical function neutral_diffusion_unit_tests() (/0.,10.,20.,30./), (/20.,16.,12.,8./), (/0.,0.,0.,0./), & ! Right positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS PiLRo, PiRLo, KoL, KoR, hEff) - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(3, KoL, KoR, PiLRo, PiRLo, hEff, & + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, KoL, KoR, PiLRo, PiRLo, hEff, & (/1,1,2,2,3,3,3,3/), & ! kL (/1,1,1,2,2,3,3,3/), & ! kR (/0.,0.5,0.,0.5,0.,0.5,1.,1./), & ! pL (/0.,0.,0.5,0.,0.5,0.,0.5,1./), & ! pR (/0.,5.,5.,5.,5.,5.,0./), & ! hEff 'Right column slightly cooler') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(8, & + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(v, 8, & absolute_positions(3, (/0.,10.,20.,30./), KoL, PiLRo), & (/0.,5.,10.,15.,20.,25.,30.,30./), '... left positions') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(8, & + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(v, 8, & absolute_positions(3, (/0.,10.,20.,30./), KoR, PiRLo), & (/0.,0.,5.,10.,15.,20.,25.,30./), '... right positions') @@ -1159,7 +1160,7 @@ logical function neutral_diffusion_unit_tests() (/0.,10.,20.,30./), (/24.,20.,16.,12./), (/0.,0.,0.,0./), & ! Right positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS PiLRo, PiRLo, KoL, KoR, hEff) - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(3, KoL, KoR, PiLRo, PiRLo, hEff, & + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, KoL, KoR, PiLRo, PiRLo, hEff, & (/1,1,1,2,2,3,3,3/), & ! kL (/1,1,2,2,3,3,3,3/), & ! kR (/0.,0.,0.5,0.,0.5,0.,0.5,1./), & ! pL @@ -1174,7 +1175,7 @@ logical function neutral_diffusion_unit_tests() (/0.,10.,20.,30./), (/16.,12.,8.,4./), (/0.,0.,0.,0./), & ! Right positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS PiLRo, PiRLo, KoL, KoR, hEff) - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(3, KoL, KoR, PiLRo, PiRLo, hEff, & + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, KoL, KoR, PiLRo, PiRLo, hEff, & (/1,2,2,3,3,3,3,3/), & ! kL (/1,1,1,1,2,2,3,3/), & ! kR (/0.,0.,0.5,0.,0.5,1.,1.,1./), & ! pL @@ -1189,7 +1190,7 @@ logical function neutral_diffusion_unit_tests() (/0.,10.,20.,30./), (/9.,7.,5.,3./), (/0.,0.,0.,0./), & ! Right positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS PiLRo, PiRLo, KoL, KoR, hEff) - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(3, KoL, KoR, PiLRo, PiRLo, hEff, & + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, KoL, KoR, PiLRo, PiRLo, hEff, & (/1,2,3,3,3,3,3,3/), & ! kL (/1,1,1,1,1,2,3,3/), & ! kR (/0.,0.,0.,1.,1.,1.,1.,1./), & ! pL @@ -1204,7 +1205,7 @@ logical function neutral_diffusion_unit_tests() (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Right positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS PiLRo, PiRLo, KoL, KoR, hEff) - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(3, KoL, KoR, PiLRo, PiRLo, hEff, & + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, KoL, KoR, PiLRo, PiRLo, hEff, & (/1,2,3,3,3,3,3,3/), & ! kL (/1,1,1,1,2,3,3,3/), & ! kR (/0.,0.,0.,0.,0.,1.,1.,1./), & ! pL @@ -1219,7 +1220,7 @@ logical function neutral_diffusion_unit_tests() (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Right positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS PiLRo, PiRLo, KoL, KoR, hEff) - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(3, KoL, KoR, PiLRo, PiRLo, hEff, & + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, KoL, KoR, PiLRo, PiRLo, hEff, & (/1,1,2,2,3,3,3,3/), & ! kL (/1,1,2,2,3,3,3,3/), & ! kR (/0.,0.,0.,0.,0.,0.,1.,1./), & ! pL @@ -1234,7 +1235,7 @@ logical function neutral_diffusion_unit_tests() (/0.,10.,20.,30./), (/10.,14.,12.,4./), (/0.,0.,0.,0./), & ! Right positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS PiLRo, PiRLo, KoL, KoR, hEff) - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(3, KoL, KoR, PiLRo, PiRLo, hEff, & + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, KoL, KoR, PiLRo, PiRLo, hEff, & (/1,2,3,3,3,3,3,3/), & ! kL (/1,1,1,2,3,3,3,3/), & ! kR (/0.,0.,0.,0.,0.,0.,.75,1./), & ! pL @@ -1249,7 +1250,7 @@ logical function neutral_diffusion_unit_tests() (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Right positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS PiLRo, PiRLo, KoL, KoR, hEff) - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(3, KoL, KoR, PiLRo, PiRLo, hEff, & + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, KoL, KoR, PiLRo, PiRLo, hEff, & (/1,1,1,2,3,3,3,3/), & ! kL (/1,2,3,3,3,3,3,3/), & ! kR (/0.,0.,0.,0.,0.,0.25,1.,1./), & ! pL @@ -1264,7 +1265,7 @@ logical function neutral_diffusion_unit_tests() (/0.,10.,20.,30./), (/10.,14.,12.,4./), (/0.,0.,0.,0./), & ! Right positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS PiLRo, PiRLo, KoL, KoR, hEff) - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(3, KoL, KoR, PiLRo, PiRLo, hEff, & + neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, KoL, KoR, PiLRo, PiRLo, hEff, & (/1,1,1,1,2,3,3,3/), & ! kL (/1,2,3,3,3,3,3,3/), & ! kR (/0.,0.,0.,0.,0.,0.,0.75,1./), & ! pL @@ -1272,236 +1273,240 @@ logical function neutral_diffusion_unit_tests() (/0.,0.,0.,0.,0.,6.,0./), & ! hEff 'Two unstable mixed layers') - write(*,'(a)') '==========================================================' - - contains - - !> Returns true if a test of fv_diff() fails, and conditionally writes results to stream - logical function test_fv_diff(hkm1, hk, hkp1, Skm1, Sk, Skp1, Ptrue, title) - real, intent(in) :: hkm1 !< Left cell width - real, intent(in) :: hk !< Center cell width - real, intent(in) :: hkp1 !< Right cell width - real, intent(in) :: Skm1 !< Left cell average value - real, intent(in) :: Sk !< Center cell average value - real, intent(in) :: Skp1 !< Right cell average value - real, intent(in) :: Ptrue !< True answer (Pa) - character(len=*), intent(in) :: title !< Title for messages - - ! Local variables - integer :: stdunit - real :: Pret - - Pret = fv_diff(hkm1, hk, hkp1, Skm1, Sk, Skp1) - test_fv_diff = (Pret /= Ptrue) - - if (test_fv_diff .or. verbosity>5) then - stdunit = 6 - if (test_fv_diff) stdunit = 0 ! In case of wrong results, write to error stream - write(stdunit,'(a)') title - if (test_fv_diff) then - write(stdunit,'(2(x,a,f20.16),x,a)') 'pRet=',Pret,'pTrue=',Ptrue,'WRONG!' - else - write(stdunit,'(2(x,a,f20.16))') 'pRet=',Pret,'pTrue=',Ptrue - endif + if (.not. neutral_diffusion_unit_tests) write(*,*) 'Pass' + +end function neutral_diffusion_unit_tests + +!> Returns true if a test of fv_diff() fails, and conditionally writes results to stream +logical function test_fv_diff(verbose, hkm1, hk, hkp1, Skm1, Sk, Skp1, Ptrue, title) + logical, intent(in) :: verbose !< If true, write results to stdout + real, intent(in) :: hkm1 !< Left cell width + real, intent(in) :: hk !< Center cell width + real, intent(in) :: hkp1 !< Right cell width + real, intent(in) :: Skm1 !< Left cell average value + real, intent(in) :: Sk !< Center cell average value + real, intent(in) :: Skp1 !< Right cell average value + real, intent(in) :: Ptrue !< True answer (Pa) + character(len=*), intent(in) :: title !< Title for messages + + ! Local variables + integer :: stdunit + real :: Pret + + Pret = fv_diff(hkm1, hk, hkp1, Skm1, Sk, Skp1) + test_fv_diff = (Pret /= Ptrue) + + if (test_fv_diff .or. verbose) then + stdunit = 6 + if (test_fv_diff) stdunit = 0 ! In case of wrong results, write to error stream + write(stdunit,'(a)') title + if (test_fv_diff) then + write(stdunit,'(2(x,a,f20.16),x,a)') 'pRet=',Pret,'pTrue=',Ptrue,'WRONG!' + else + write(stdunit,'(2(x,a,f20.16))') 'pRet=',Pret,'pTrue=',Ptrue endif + endif - end function test_fv_diff - - !> Returns true if a test of fvlsq_slope() fails, and conditionally writes results to stream - logical function test_fvlsq_slope(hkm1, hk, hkp1, Skm1, Sk, Skp1, Ptrue, title) - real, intent(in) :: hkm1 !< Left cell width - real, intent(in) :: hk !< Center cell width - real, intent(in) :: hkp1 !< Right cell width - real, intent(in) :: Skm1 !< Left cell average value - real, intent(in) :: Sk !< Center cell average value - real, intent(in) :: Skp1 !< Right cell average value - real, intent(in) :: Ptrue !< True answer (Pa) - character(len=*), intent(in) :: title !< Title for messages - - ! Local variables - integer :: stdunit - real :: Pret - - Pret = fvlsq_slope(hkm1, hk, hkp1, Skm1, Sk, Skp1) - test_fvlsq_slope = (Pret /= Ptrue) - - if (test_fvlsq_slope .or. verbosity>5) then - stdunit = 6 - if (test_fvlsq_slope) stdunit = 0 ! In case of wrong results, write to error stream - write(stdunit,'(a)') title - if (test_fvlsq_slope) then - write(stdunit,'(2(x,a,f20.16),x,a)') 'pRet=',Pret,'pTrue=',Ptrue,'WRONG!' - else - write(stdunit,'(2(x,a,f20.16))') 'pRet=',Pret,'pTrue=',Ptrue - endif +end function test_fv_diff + +!> Returns true if a test of fvlsq_slope() fails, and conditionally writes results to stream +logical function test_fvlsq_slope(verbose, hkm1, hk, hkp1, Skm1, Sk, Skp1, Ptrue, title) + logical, intent(in) :: verbose !< If true, write results to stdout + real, intent(in) :: hkm1 !< Left cell width + real, intent(in) :: hk !< Center cell width + real, intent(in) :: hkp1 !< Right cell width + real, intent(in) :: Skm1 !< Left cell average value + real, intent(in) :: Sk !< Center cell average value + real, intent(in) :: Skp1 !< Right cell average value + real, intent(in) :: Ptrue !< True answer (Pa) + character(len=*), intent(in) :: title !< Title for messages + + ! Local variables + integer :: stdunit + real :: Pret + + Pret = fvlsq_slope(hkm1, hk, hkp1, Skm1, Sk, Skp1) + test_fvlsq_slope = (Pret /= Ptrue) + + if (test_fvlsq_slope .or. verbose) then + stdunit = 6 + if (test_fvlsq_slope) stdunit = 0 ! In case of wrong results, write to error stream + write(stdunit,'(a)') title + if (test_fvlsq_slope) then + write(stdunit,'(2(x,a,f20.16),x,a)') 'pRet=',Pret,'pTrue=',Ptrue,'WRONG!' + else + write(stdunit,'(2(x,a,f20.16))') 'pRet=',Pret,'pTrue=',Ptrue endif + endif - end function test_fvlsq_slope - - !> Returns true if a test of interpolate_for_nondim_position() fails, and conditionally writes results to stream - logical function test_ifndp(rhoNeg, Pneg, rhoPos, Ppos, Ptrue, title) - real, intent(in) :: rhoNeg !< Lighter density (kg/m3) - real, intent(in) :: Pneg !< Interface position of lighter density (Pa) - real, intent(in) :: rhoPos !< Heavier density (kg/m3) - real, intent(in) :: Ppos !< Interface position of heavier density (Pa) - real, intent(in) :: Ptrue !< True answer (Pa) - character(len=*), intent(in) :: title !< Title for messages - - ! Local variables - integer :: stdunit - real :: Pret - - Pret = interpolate_for_nondim_position(rhoNeg, Pneg, rhoPos, Ppos) - test_ifndp = (Pret /= Ptrue) - - if (test_ifndp .or. verbosity>5) then - stdunit = 6 - if (test_ifndp) stdunit = 0 ! In case of wrong results, write to error stream - write(stdunit,'(a)') title - if (test_ifndp) then - write(stdunit,'(4(x,a,f20.16),2(x,a,1pe22.15),x,a)') 'r1=',rhoNeg,'p1=',Pneg,'r2=',rhoPos,'p2=',Ppos,'pRet=',Pret,'pTrue=',Ptrue,'WRONG!' - else - write(stdunit,'(4(x,a,f20.16),2(x,a,1pe22.15))') 'r1=',rhoNeg,'p1=',Pneg,'r2=',rhoPos,'p2=',Ppos,'pRet=',Pret,'pTrue=',Ptrue - endif +end function test_fvlsq_slope + +!> Returns true if a test of interpolate_for_nondim_position() fails, and conditionally writes results to stream +logical function test_ifndp(verbose, rhoNeg, Pneg, rhoPos, Ppos, Ptrue, title) + logical, intent(in) :: verbose !< If true, write results to stdout + real, intent(in) :: rhoNeg !< Lighter density (kg/m3) + real, intent(in) :: Pneg !< Interface position of lighter density (Pa) + real, intent(in) :: rhoPos !< Heavier density (kg/m3) + real, intent(in) :: Ppos !< Interface position of heavier density (Pa) + real, intent(in) :: Ptrue !< True answer (Pa) + character(len=*), intent(in) :: title !< Title for messages + + ! Local variables + integer :: stdunit + real :: Pret + + Pret = interpolate_for_nondim_position(rhoNeg, Pneg, rhoPos, Ppos) + test_ifndp = (Pret /= Ptrue) + + if (test_ifndp .or. verbose) then + stdunit = 6 + if (test_ifndp) stdunit = 0 ! In case of wrong results, write to error stream + write(stdunit,'(a)') title + if (test_ifndp) then + write(stdunit,'(4(x,a,f20.16),2(x,a,1pe22.15),x,a)') 'r1=',rhoNeg,'p1=',Pneg,'r2=',rhoPos,'p2=',Ppos,'pRet=',Pret,'pTrue=',Ptrue,'WRONG!' + else + write(stdunit,'(4(x,a,f20.16),2(x,a,1pe22.15))') 'r1=',rhoNeg,'p1=',Pneg,'r2=',rhoPos,'p2=',Ppos,'pRet=',Pret,'pTrue=',Ptrue endif + endif - end function test_ifndp +end function test_ifndp - !> Returns true if comparison of Po and Ptrue fails, and conditionally writes results to stream - logical function test_data1d(nk, Po, Ptrue, title) - integer, intent(in) :: nk !< Number of layers - real, dimension(nk), intent(in) :: Po !< Calculated answer - real, dimension(nk), intent(in) :: Ptrue !< True answer - character(len=*), intent(in) :: title !< Title for messages +!> Returns true if comparison of Po and Ptrue fails, and conditionally writes results to stream +logical function test_data1d(verbose, nk, Po, Ptrue, title) + logical, intent(in) :: verbose !< If true, write results to stdout + integer, intent(in) :: nk !< Number of layers + real, dimension(nk), intent(in) :: Po !< Calculated answer + real, dimension(nk), intent(in) :: Ptrue !< True answer + character(len=*), intent(in) :: title !< Title for messages - ! Local variables - integer :: k, stdunit + ! Local variables + integer :: k, stdunit - test_data1d = .false. + test_data1d = .false. + do k = 1,nk + if (Po(k) /= Ptrue(k)) test_data1d = .true. + enddo + + if (test_data1d .or. verbose) then + stdunit = 6 + if (test_data1d) stdunit = 0 ! In case of wrong results, write to error stream + write(stdunit,'(a)') title do k = 1,nk - if (Po(k) /= Ptrue(k)) test_data1d = .true. + if (Po(k) /= Ptrue(k)) then + test_data1d = .true. + write(stdunit,'(a,i2,2(x,a,f20.16),x,a,1pe22.15,x,a)') 'k=',k,'Po=',Po(k),'Ptrue=',Ptrue(k),'err=',Po(k)-Ptrue(k),'WRONG!' + else + if (verbose) & + write(stdunit,'(a,i2,2(x,a,f20.16),x,a,1pe22.15)') 'k=',k,'Po=',Po(k),'Ptrue=',Ptrue(k),'err=',Po(k)-Ptrue(k) + endif enddo + endif - if (test_data1d .or. verbosity>5) then - stdunit = 6 - if (test_data1d) stdunit = 0 ! In case of wrong results, write to error stream - write(stdunit,'(a)') title - do k = 1,nk - if (Po(k) /= Ptrue(k)) then - test_data1d = .true. - write(stdunit,'(a,i2,2(x,a,f20.16),x,a,1pe22.15,x,a)') 'k=',k,'Po=',Po(k),'Ptrue=',Ptrue(k),'err=',Po(k)-Ptrue(k),'WRONG!' - else - if (verbosity>5) & - write(stdunit,'(a,i2,2(x,a,f20.16),x,a,1pe22.15)') 'k=',k,'Po=',Po(k),'Ptrue=',Ptrue(k),'err=',Po(k)-Ptrue(k) - endif - enddo - endif +end function test_data1d - end function test_data1d +!> Returns true if comparison of Po and Ptrue fails, and conditionally writes results to stream +logical function test_data1di(verbose, nk, Po, Ptrue, title) + logical, intent(in) :: verbose !< If true, write results to stdout + integer, intent(in) :: nk !< Number of layers + integer, dimension(nk), intent(in) :: Po !< Calculated answer + integer, dimension(nk), intent(in) :: Ptrue !< True answer + character(len=*), intent(in) :: title !< Title for messages - !> Returns true if comparison of Po and Ptrue fails, and conditionally writes results to stream - logical function test_data1di(nk, Po, Ptrue, title) - integer, intent(in) :: nk !< Number of layers - integer, dimension(nk), intent(in) :: Po !< Calculated answer - integer, dimension(nk), intent(in) :: Ptrue !< True answer - character(len=*), intent(in) :: title !< Title for messages + ! Local variables + integer :: k, stdunit - ! Local variables - integer :: k, stdunit + test_data1di = .false. + do k = 1,nk + if (Po(k) /= Ptrue(k)) test_data1di = .true. + enddo - test_data1di = .false. + if (test_data1di .or. verbose) then + stdunit = 6 + if (test_data1di) stdunit = 0 ! In case of wrong results, write to error stream + write(stdunit,'(a)') title do k = 1,nk - if (Po(k) /= Ptrue(k)) test_data1di = .true. + if (Po(k) /= Ptrue(k)) then + test_data1di = .true. + write(stdunit,'(a,i2,2(x,a,i5),x,a)') 'k=',k,'Io=',Po(k),'Itrue=',Ptrue(k),'WRONG!' + else + if (verbose) & + write(stdunit,'(a,i2,2(x,a,i5))') 'k=',k,'Io=',Po(k),'Itrue=',Ptrue(k) + endif enddo + endif - if (test_data1di .or. verbosity>5) then - stdunit = 6 - if (test_data1di) stdunit = 0 ! In case of wrong results, write to error stream - write(stdunit,'(a)') title - do k = 1,nk - if (Po(k) /= Ptrue(k)) then - test_data1di = .true. - write(stdunit,'(a,i2,2(x,a,i5),x,a)') 'k=',k,'Io=',Po(k),'Itrue=',Ptrue(k),'WRONG!' - else - if (verbosity>5) & - write(stdunit,'(a,i2,2(x,a,i5))') 'k=',k,'Io=',Po(k),'Itrue=',Ptrue(k) - endif - enddo +end function test_data1di + +!> Returns true if output of find_neutral_surface_positions() does not match correct values, and conditionally writes results to stream +logical function test_nsp(verbose, nk, KoL, KoR, pL, pR, hEff, KoL0, KoR0, pL0, pR0, hEff0, title) + logical, intent(in) :: verbose !< If true, write results to stdout + integer, intent(in) :: nk !< Number of layers + integer, dimension(2*nk+2), intent(in) :: KoL !< Index of first left interface above neutral surface + integer, dimension(2*nk+2), intent(in) :: KoR !< Index of first right interface above neutral surface + real, dimension(2*nk+2), intent(in) :: pL !< Fractional position of neutral surface within layer KoL of left column + real, dimension(2*nk+2), intent(in) :: pR !< Fractional position of neutral surface within layer KoR of right column + real, dimension(2*nk+1), intent(in) :: hEff !< Effective thickness between two neutral surfaces (Pa) + integer, dimension(2*nk+2), intent(in) :: KoL0 !< Correct value for KoL + integer, dimension(2*nk+2), intent(in) :: KoR0 !< Correct value for KoR + real, dimension(2*nk+2), intent(in) :: pL0 !< Correct value for pL + real, dimension(2*nk+2), intent(in) :: pR0 !< Correct value for pR + real, dimension(2*nk+1), intent(in) :: hEff0 !< Correct value for hEff + character(len=*), intent(in) :: title !< Title for messages + + ! Local variables + integer :: k, stdunit + logical :: this_row_failed + + test_nsp = .false. + do k = 1,2*nk+2 + test_nsp = test_nsp .or. compare_nsp_row(KoL(k), KoR(k), pL(k), pR(k), KoL0(k), KoR0(k), pL0(k), pR0(k)) + if (k < 2*nk+2) then + if (hEff(k) /= hEff0(k)) test_nsp = .true. endif + enddo - end function test_data1di - - !> Returns true if output of find_neutral_surface_positions() does not match correct values, and conditionally writes results to stream - logical function test_nsp(nk, KoL, KoR, pL, pR, hEff, KoL0, KoR0, pL0, pR0, hEff0, title) - integer, intent(in) :: nk !< Number of layers - integer, dimension(2*nk+2), intent(in) :: KoL !< Index of first left interface above neutral surface - integer, dimension(2*nk+2), intent(in) :: KoR !< Index of first right interface above neutral surface - real, dimension(2*nk+2), intent(in) :: pL !< Fractional position of neutral surface within layer KoL of left column - real, dimension(2*nk+2), intent(in) :: pR !< Fractional position of neutral surface within layer KoR of right column - real, dimension(2*nk+1), intent(in) :: hEff !< Effective thickness between two neutral surfaces (Pa) - integer, dimension(2*nk+2), intent(in) :: KoL0 !< Correct value for KoL - integer, dimension(2*nk+2), intent(in) :: KoR0 !< Correct value for KoR - real, dimension(2*nk+2), intent(in) :: pL0 !< Correct value for pL - real, dimension(2*nk+2), intent(in) :: pR0 !< Correct value for pR - real, dimension(2*nk+1), intent(in) :: hEff0 !< Correct value for hEff - character(len=*), intent(in) :: title !< Title for messages - - ! Local variables - integer :: k, stdunit - logical :: this_row_failed - - test_nsp = .false. + if (test_nsp .or. verbose) then + stdunit = 6 + if (test_nsp) stdunit = 0 ! In case of wrong results, write to error stream + write(stdunit,'(a)') title do k = 1,2*nk+2 - test_nsp = test_nsp .or. compare_nsp_row(KoL(k), KoR(k), pL(k), pR(k), KoL0(k), KoR0(k), pL0(k), pR0(k)) - if (k < 2*nk+2) then - if (hEff(k) /= hEff0(k)) test_nsp = .true. + this_row_failed = compare_nsp_row(KoL(k), KoR(k), pL(k), pR(k), KoL0(k), KoR0(k), pL0(k), pR0(k)) + if (this_row_failed) then + write(stdunit,10) k,KoL(k),pL(k),KoR(k),pR(k),' <-- WRONG!' + write(stdunit,10) k,KoL0(k),pL0(k),KoR0(k),pR0(k),' <-- should be this' + else + write(stdunit,10) k,KoL(k),pL(k),KoR(k),pR(k) endif - enddo - - if (test_nsp .or. verbosity>5) then - stdunit = 6 - if (test_nsp) stdunit = 0 ! In case of wrong results, write to error stream - write(stdunit,'(a)') title - do k = 1,2*nk+2 - this_row_failed = compare_nsp_row(KoL(k), KoR(k), pL(k), pR(k), KoL0(k), KoR0(k), pL0(k), pR0(k)) - if (this_row_failed) then - write(stdunit,10) k,KoL(k),pL(k),KoR(k),pR(k),' <-- WRONG!' - write(stdunit,10) k,KoL0(k),pL0(k),KoR0(k),pR0(k),' <-- should be this' + if (k < 2*nk+2) then + if (hEff(k) /= hEff0(k)) then + write(stdunit,'(i3,8x,"layer hEff =",2(f20.16,a))') k,hEff(k)," .neq. ",hEff0(k),' <-- WRONG!' else - write(stdunit,10) k,KoL(k),pL(k),KoR(k),pR(k) - endif - if (k < 2*nk+2) then - if (hEff(k) /= hEff0(k)) then - write(stdunit,'(i3,8x,"layer hEff =",2(f20.16,a))') k,hEff(k)," .neq. ",hEff0(k),' <-- WRONG!' - else - write(stdunit,'(i3,8x,"layer hEff =",f20.16)') k,hEff(k) - endif + write(stdunit,'(i3,8x,"layer hEff =",f20.16)') k,hEff(k) endif - enddo - endif - -10 format("ks=",i3," kL=",i3," pL=",f20.16," kR=",i3," pR=",f20.16,a) - end function test_nsp - - !> Compares a single row, k, of output from find_neutral_surface_positions() - logical function compare_nsp_row(KoL, KoR, pL, pR, KoL0, KoR0, pL0, pR0) - integer, intent(in) :: KoL !< Index of first left interface above neutral surface - integer, intent(in) :: KoR !< Index of first right interface above neutral surface - real, intent(in) :: pL !< Fractional position of neutral surface within layer KoL of left column - real, intent(in) :: pR !< Fractional position of neutral surface within layer KoR of right column - integer, intent(in) :: KoL0 !< Correct value for KoL - integer, intent(in) :: KoR0 !< Correct value for KoR - real, intent(in) :: pL0 !< Correct value for pL - real, intent(in) :: pR0 !< Correct value for pR - - compare_nsp_row = .false. - if (KoL /= KoL0) compare_nsp_row = .true. - if (KoR /= KoR0) compare_nsp_row = .true. - if (pL /= pL0) compare_nsp_row = .true. - if (pR /= pR0) compare_nsp_row = .true. - end function compare_nsp_row + endif + enddo + endif -end function neutral_diffusion_unit_tests +10 format("ks=",i3," kL=",i3," pL=",f20.16," kR=",i3," pR=",f20.16,a) +end function test_nsp + +!> Compares a single row, k, of output from find_neutral_surface_positions() +logical function compare_nsp_row(KoL, KoR, pL, pR, KoL0, KoR0, pL0, pR0) + integer, intent(in) :: KoL !< Index of first left interface above neutral surface + integer, intent(in) :: KoR !< Index of first right interface above neutral surface + real, intent(in) :: pL !< Fractional position of neutral surface within layer KoL of left column + real, intent(in) :: pR !< Fractional position of neutral surface within layer KoR of right column + integer, intent(in) :: KoL0 !< Correct value for KoL + integer, intent(in) :: KoR0 !< Correct value for KoR + real, intent(in) :: pL0 !< Correct value for pL + real, intent(in) :: pR0 !< Correct value for pR + + compare_nsp_row = .false. + if (KoL /= KoL0) compare_nsp_row = .true. + if (KoR /= KoR0) compare_nsp_row = .true. + if (pL /= pL0) compare_nsp_row = .true. + if (pR /= pR0) compare_nsp_row = .true. +end function compare_nsp_row !> Deallocates neutral_diffusion control structure subroutine neutral_diffusion_end(CS) diff --git a/src/tracer/ideal_age_example.F90 b/src/tracer/ideal_age_example.F90 index dd4200d61b..7e544ec48f 100644 --- a/src/tracer/ideal_age_example.F90 +++ b/src/tracer/ideal_age_example.F90 @@ -313,6 +313,7 @@ subroutine initialize_ideal_age_tracer(restart, day, G, GV, h, diag, OBC, CS, & character(len=48) :: units ! The dimensions of the variable. character(len=48) :: flux_units ! The units for age tracer fluxes, either ! years m3 s-1 or years kg s-1. + character(len=72) :: cmorname ! The CMOR name of that variable. logical :: OK integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz, m integer :: IsdB, IedB, JsdB, JedB @@ -379,11 +380,15 @@ subroutine initialize_ideal_age_tracer(restart, day, G, GV, h, diag, OBC, CS, & else ; flux_units = "years kg s-1" ; endif do m=1,CS%ntr - ! Register the tracer for the restart file. call query_vardesc(CS%tr_desc(m), name, units=units, longname=longname, & - caller="initialize_ideal_age_tracer") - CS%id_tracer(m) = register_diag_field("ocean_model", trim(name), CS%diag%axesTL, & + cmor_field_name=cmorname, caller="initialize_ideal_age_tracer") + if (len_trim(cmorname)==0) then + CS%id_tracer(m) = register_diag_field("ocean_model", trim(name), CS%diag%axesTL, & day, trim(longname) , trim(units)) + else + CS%id_tracer(m) = register_diag_field("ocean_model", trim(name), CS%diag%axesTL, & + day, trim(longname) , trim(units), cmor_field_name=cmorname) + endif CS%id_tr_adx(m) = register_diag_field("ocean_model", trim(name)//"_adx", & CS%diag%axesCuL, day, trim(longname)//" advective zonal flux" , & trim(flux_units))