diff --git a/.codecov.yml b/.codecov.yml index 576633bf6a..05fe474ab3 100644 --- a/.codecov.yml +++ b/.codecov.yml @@ -6,3 +6,6 @@ coverage: patch: default: threshold: 100% +comment: + # This must be set to the number of test cases (TCs) + after_n_builds: 8 diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 39b63c8f85..5a05694fef 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -32,13 +32,11 @@ setup: - git clone --recursive http://gitlab.gfdl.noaa.gov/ogrp/Gaea-stats-MOM6-examples.git tests && cd tests # Install / update testing scripts - git clone https://github.com/adcroft/MRS.git MRS - - (cd MRS ; git checkout xanadu-fms) # Update MOM6-examples and submodules - (cd MOM6-examples && git checkout . && git checkout dev/gfdl && git pull && git submodule init && git submodule update) - (cd MOM6-examples/src/MOM6 && git submodule update) - test -d MOM6-examples/src/LM3 || make -f MRS/Makefile.clone clone_gfdl -s - make -f MRS/Makefile.clone MOM6-examples/.datasets -s - #- (cd MOM6-examples/src/mkmf && git pull https://github.com/adcroft/mkmf.git add_coverage_mode) - env > gitlab_session.log # Cache everything under tests to unpack for each subsequent stage - cd ../ ; time tar zcf $CACHE_DIR/tests_$CI_PIPELINE_ID.tgz tests diff --git a/.testing/Makefile b/.testing/Makefile index 99672268c3..d38189667b 100644 --- a/.testing/Makefile +++ b/.testing/Makefile @@ -119,28 +119,28 @@ build/%/MOM6: build/%/Makefile $(FMS)/lib/libfms.a build/%/Makefile: build/%/path_names cp $(MKMF_TEMPLATE) $(@D) cd $(@D) && $(MKMF) \ - -t $(notdir $(MKMF_TEMPLATE)) \ - -o '-I ../../$(DEPS)/fms/build' \ - -p MOM6 \ - -l '../../$(DEPS)/fms/lib/libfms.a' \ - -c $(MKMF_CPP) \ - path_names + -t $(notdir $(MKMF_TEMPLATE)) \ + -o '-I ../../$(DEPS)/fms/build' \ + -p MOM6 \ + -l '../../$(DEPS)/fms/lib/libfms.a' \ + -c $(MKMF_CPP) \ + path_names # NOTE: These path_names rules could be merged build/target/path_names: $(LIST_PATHS) $(TARGET_CODEBASE) $(TARGET_SOURCE) mkdir -p $(@D) cd $(@D) && $(LIST_PATHS) -l \ - ../../$(TARGET_CODEBASE)/src \ - ../../$(TARGET_CODEBASE)/config_src/solo_driver \ - ../../$(TARGET_CODEBASE)/$(GRID_SRC) + ../../$(TARGET_CODEBASE)/src \ + ../../$(TARGET_CODEBASE)/config_src/solo_driver \ + ../../$(TARGET_CODEBASE)/$(GRID_SRC) build/%/path_names: $(LIST_PATHS) $(MOM_SOURCE) mkdir -p $(@D) cd $(@D) && $(LIST_PATHS) -l \ - ../../../src \ - ../../../config_src/solo_driver \ - ../../../$(GRID_SRC) + ../../../src \ + ../../../config_src/solo_driver \ + ../../../$(GRID_SRC) # Target repository for regression tests $(TARGET_CODEBASE): @@ -158,10 +158,10 @@ $(FMS)/lib/libfms.a: $(FMS)/build/Makefile $(FMS)/build/Makefile: $(FMS)/build/path_names cp $(MKMF_TEMPLATE) $(@D) cd $(@D) && $(MKMF) \ - -t $(notdir $(MKMF_TEMPLATE)) \ - -p ../lib/libfms.a \ - -c $(MKMF_CPP) \ - path_names + -t $(notdir $(MKMF_TEMPLATE)) \ + -p ../lib/libfms.a \ + -c $(MKMF_CPP) \ + path_names $(FMS)/build/path_names: $(LIST_PATHS) $(FMS)/src $(FMS_SOURCE) mkdir -p $(@D) @@ -202,18 +202,38 @@ test.repros: $(foreach c,$(CONFIGS),$(c).repro $(c).repro.diag) test.openmps: $(foreach c,$(CONFIGS),$(c).openmp $(c).openmp.diag) test.nans: $(foreach c,$(CONFIGS),$(c).nan $(c).nan.diag) test.dims: $(foreach c,$(CONFIGS),$(foreach d,$(DIMS),$(c).dim.$(d) $(c).dim.$(d).diag)) - test.regressions: $(foreach c,$(CONFIGS),$(c).regression $(c).regression.diag) - ! ls -1 results/*/*.reg -define CMP_RULE -.PRECIOUS: $(foreach b,$(2),results/%/ocean.stats.$(b)) -%.$(1): $(foreach b,$(2),results/%/ocean.stats.$(b)) - cmp $$^ || diff $$^ +# Color highlights for test results +RED=\033[0;31m +GREEN=\033[0;32m +RESET=\033[0m + +DONE=${GREEN}DONE${RESET} +PASS=${GREEN}PASS${RESET} +FAIL=${RED}FAIL${RESET} -.PRECIOUS: $(foreach b,$(2),results/%/chksum_diag.$(b)) -%.$(1).diag: $(foreach b,$(2),results/%/chksum_diag.$(b)) - cmp $$^ || diff $$^ +# Comparison rules +# $(1): Test type (grid, layout, &c.) +# $(2): Comparison targets (symmetric asymmetric, symmetric layout, &c.) +define CMP_RULE +.PRECIOUS: $(foreach b,$(2),work/%/$(b)/ocean.stats) +%.$(1): $(foreach b,$(2),work/%/$(b)/ocean.stats) + @cmp $$^ || !( \ + mkdir -p results/$$*; \ + (diff $$^ | tee results/$$*/ocean.stats.$(1).diff | head) ; \ + echo -e "${FAIL}: Solutions $$*.$(1) have changed." \ + ) + @echo -e "${PASS}: Solutions $$*.$(1) agree." + +.PRECIOUS: $(foreach b,$(2),work/%/$(b)/chksum_diag) +%.$(1).diag: $(foreach b,$(2),work/%/$(b)/chksum_diag) + @cmp $$^ || !( \ + mkdir -p results/$$*; \ + (diff $$^ | tee results/$$*/chksum_diag.$(1).diff | head) ; \ + echo -e "${FAIL}: Diagnostics $$*.$(1).diag have changed." \ + ) + @echo -e "${PASS}: Diagnostics $$*.$(1).diag agree." endef $(eval $(call CMP_RULE,grid,symmetric asymmetric)) @@ -223,29 +243,31 @@ $(eval $(call CMP_RULE,repro,symmetric repro)) $(eval $(call CMP_RULE,openmp,symmetric openmp)) $(eval $(call CMP_RULE,nan,symmetric nan)) $(foreach d,$(DIMS),$(eval $(call CMP_RULE,dim.$(d),symmetric dim.$(d)))) +$(eval $(call CMP_RULE,regression,symmetric target)) # Custom comparison rules -.PRECIOUS: $(foreach b,symmetric restart target,results/%/ocean.stats.$(b)) - # Restart tests only compare the final stat record -%.restart: $(foreach b,symmetric restart,results/%/ocean.stats.$(b)) - cmp $(foreach f,$^,<(tr -s ' ' < $(f) | cut -d ' ' -f3- | tail -n 1)) \ - || diff $^ +.PRECIOUS: $(foreach b,symmetric restart target,work/%/$(b)/ocean.stats) +%.restart: $(foreach b,symmetric restart,work/%/$(b)/ocean.stats) + #cmp $(foreach f,$^,<(tr -s ' ' < $(f) | cut -d ' ' -f3- | tail -n 1)) \ + # || diff $^ + @cmp $(foreach f,$^,<(tr -s ' ' < $(f) | cut -d ' ' -f3- | tail -n 1)) \ + || !( \ + mkdir -p results/$*; \ + (diff $$^ | tee results/$*/chksum_diag.restart.diff | head) ; \ + echo -e "${FAIL}: Diagnostics $*.restart.diag have changed." \ + ) + @echo -e "${PASS}: Diagnostics $*.restart.diag agree." # TODO: chksum_diag parsing of restart files -# All regression tests must be completed when considering answer changes -%.regression: $(foreach b,symmetric target,results/%/ocean.stats.$(b)) - cmp $^ || (diff $^ > $<.reg || true) - -%.regression.diag: $(foreach b,symmetric target,results/%/chksum_diag.$(b)) - cmp $^ || (diff $^ > $<.reg || true) #--- # Test run output files # Generalized MPI environment variable support +# XXX: Using `-env` in the MPICH test can erroneously producing an `nv` file. # $(1): Environment variables ifeq ($(shell $(MPIRUN) -x tmp=1 true 2> /dev/null ; echo $$?), 0) MPIRUN_CMD=$(MPIRUN) $(if $(1),-x $(1),) @@ -255,7 +277,8 @@ else MPIRUN_CMD=$(1) $(MPIRUN) endif -# Rule to build results//{ocean.stats,chksum_diag}. + +# Rule to build work//{ocean.stats,chksum_diag}. # $(1): Test configuration name # $(2): Executable type # $(3): Enable coverage flag @@ -263,25 +286,33 @@ endif # $(5): Environment variables # $(6): Number of MPI ranks define STAT_RULE -results/%/ocean.stats.$(1): build/$(2)/MOM6 +work/%/$(1)/ocean.stats work/%/$(1)/chksum_diag: build/$(2)/MOM6 + @echo "Running test $$*.$(1)..." if [ $(3) ]; then find build/$(2) -name *.gcda -exec rm -f '{}' \; ; fi - mkdir -p work/$$*/$(1) - cp -rL $$*/* work/$$*/$(1) - cd work/$$*/$(1) && if [ -f Makefile ]; then make; fi - mkdir -p work/$$*/$(1)/RESTART - echo -e "$(4)" > work/$$*/$(1)/MOM_override - cd work/$$*/$(1) && $$(call MPIRUN_CMD,$(5)) -n $(6) ../../../$$< 2> debug.out > std.out \ - || ! sed 's/^/$$*.$(1): /' std.out debug.out \ - && sed 's/^/$$*.$(1): /' std.out - mkdir -p $$(@D) - cp work/$$*/$(1)/ocean.stats $$@ - if [ $(3) ]; then cd .. && bash <(curl -s https://codecov.io/bash) -n $$@; fi - -results/%/chksum_diag.$(1): results/%/ocean.stats.$(1) mkdir -p $$(@D) - cp work/$$*/$(1)/chksum_diag $$@ + cp -rL $$*/* $$(@D) + cd $$(@D) && if [ -f Makefile ]; then make; fi + mkdir -p $$(@D)/RESTART + echo -e "$(4)" > $$(@D)/MOM_override + cd $$(@D) \ + && $$(call MPIRUN_CMD,$(5)) -n $(6) ../../../$$< 2> std.err > std.out \ + || !( \ + mkdir -p ../../../results/$$*/ ; \ + cat std.out | tee ../../../results/$$*/std.$(1).out | tail ; \ + cat std.err | tee ../../../results/$$*/std.$(1).err | tail ; \ + rm ocean.stats chksum_diag ; \ + echo -e "${FAIL}: $$*.$(1) failed at runtime." \ + ) + @echo -e "${DONE}: $$*.$(1); no runtime errors." + if [ $(3) ]; then \ + mkdir -p results/$$* ; \ + bash <(curl -s https://codecov.io/bash) -n $$@ \ + > work/$$*/codecov.$(1).out \ + 2> work/$$*/codecov.$(1).err ; \ + fi endef + # Define $(,) as comma escape character , := , @@ -300,50 +331,80 @@ $(eval $(call STAT_RULE,dim.z,symmetric,,Z_RESCALE_POWER=11,,1)) $(eval $(call STAT_RULE,dim.q,symmetric,,Q_RESCALE_POWER=11,,1)) $(eval $(call STAT_RULE,dim.r,symmetric,,R_RESCALE_POWER=11,,1)) + # Restart tests require significant preprocessing, and are handled separately. -results/%/ocean.stats.restart: build/symmetric/MOM6 - rm -rf work/$*/restart - mkdir -p work/$*/restart - cp -rL $*/* work/$*/restart +work/%/restart/ocean.stats: build/symmetric/MOM6 + rm -rf $(@D) + mkdir -p $(@D) + cp -rL $*/* $(@D) cd work/$*/restart && if [ -f Makefile ]; then make; fi - mkdir -p work/$*/restart/RESTART + mkdir -p $(@D)/RESTART # Generate the half-period input namelist - # TODO: Assumes runtime set by DAYMAX, will fail if set by input.nml - cd work/$*/restart \ - && daymax=$$(grep DAYMAX MOM_input | cut -d '!' -f 1 | cut -d '=' -f 2 | xargs) \ - && timeunit=$$(grep TIMEUNIT MOM_input | cut -d '!' -f 1 | cut -d '=' -f 2 | xargs) \ - && if [ -z "$${timeunit}" ]; then timeunit="8.64e4"; fi \ - && printf -v timeunit_int "%.f" "$${timeunit}" \ - && halfperiod=$$(printf "%.f" $$(bc <<< "scale=10; 0.5 * $${daymax} * $${timeunit_int}")) \ - && printf "\n&ocean_solo_nml\n seconds = $${halfperiod}\n/\n" >> input.nml + # TODO: Assumes that runtime set by DAYMAX, will fail if set by input.nml + cd $(@D) \ + && daymax=$$(grep DAYMAX MOM_input | cut -d '!' -f 1 | cut -d '=' -f 2 | xargs) \ + && timeunit=$$(grep TIMEUNIT MOM_input | cut -d '!' -f 1 | cut -d '=' -f 2 | xargs) \ + && if [ -z "$${timeunit}" ]; then timeunit="8.64e4"; fi \ + && printf -v timeunit_int "%.f" "$${timeunit}" \ + && halfperiod=$$(printf "%.f" $$(bc <<< "scale=10; 0.5 * $${daymax} * $${timeunit_int}")) \ + && printf "\n&ocean_solo_nml\n seconds = $${halfperiod}\n/\n" >> input.nml # Run the first half-period - cd work/$*/restart && $(MPIRUN) -n 1 ../../../$< 2> debug1.out > std1.out \ - || ! sed 's/^/$*.restart1: /' std1.out debug1.out \ - && sed 's/^/$*.restart1: /' std1.out + cd $(@D) && $(MPIRUN) -n 1 ../../../$< 2> std1.err > std1.out \ + || !( \ + cat std1.out | tee ../../../results/$*/std.restart1.out | tail ; \ + cat std1.err | tee ../../../results/$*/std.restart1.err | tail ; \ + echo -e "${FAIL}: $*.restart failed at runtime." \ + ) # Setup the next inputs - cd work/$*/restart && rm -rf INPUT && mv RESTART INPUT - mkdir work/$*/restart/RESTART - cd work/$*/restart && sed -i -e "s/input_filename *= *'n'/input_filename = 'r'/g" input.nml + cd $(@D) && rm -rf INPUT && mv RESTART INPUT + mkdir $(@D)/RESTART + cd $(@D) && sed -i -e "s/input_filename *= *'n'/input_filename = 'r'/g" input.nml # Run the second half-period - cd work/$*/restart && $(MPIRUN) -n 1 ../../../$< 2> debug2.out > std2.out \ - || ! sed 's/^/$*.restart2: /' std2.out debug2.out \ - && sed 's/^/$*.restart2: /' std2.out - # Archive the results and cleanup - mkdir -p $(@D) - cp work/$*/restart/ocean.stats $@ + cd $(@D) && $(MPIRUN) -n 1 ../../../$< 2> std2.err > std2.out \ + || !( \ + cat std2.out | tee ../../../results/$*/std.restart2.out | tail ; \ + cat std2.err | tee ../../../results/$*/std.restart2.err | tail ; \ + echo -e "${FAIL}: $*.restart failed at runtime." \ + ) # TODO: Restart checksum diagnostics +#--- +# Not a true rule; only call this after `make test` to summarize test results. +.PHONY: test.summary +test.summary: + @if ls results/*/* &> /dev/null; then \ + if ls results/*/std.*.err &> /dev/null; then \ + echo "The following tests failed to complete:" ; \ + ls results/*/std.*.out \ + | awk '{split($$0,a,"/"); split(a[3],t,"."); v=t[2]; if(length(t)>3) v=v"."t[3]; print a[2],":",v}'; \ + fi; \ + if ls results/*/ocean.stats.*.diff &> /dev/null; then \ + echo "The following tests report solution regressions:" ; \ + ls results/*/ocean.stats.*.diff \ + | awk '{split($$0,a,"/"); split(a[3],t,"."); v=t[3]; if(length(t)>4) v=v"."t[4]; print a[2],":",v}'; \ + fi; \ + if ls results/*/chksum_diag.*.diff &> /dev/null; then \ + echo "The following tests report diagnostic regressions:" ; \ + ls results/*/chksum_diag.*.diff \ + | awk '{split($$0,a,"/"); split(a[3],t,"."); v=t[2]; if(length(t)>3) v=v"."t[3]; print a[2],":",v}'; \ + fi; \ + false ; \ + else \ + echo -e "${PASS}: All tests passed!"; \ + fi + + #---- +# NOTE: These tests assert that we are in the .testing directory. + .PHONY: clean clean: clean.stats - @# Assert that we are in .testing for recursive delete @[ $$(basename $$(pwd)) = .testing ] rm -rf build .PHONY: clean.stats clean.stats: - @# Assert that we are in .testing for recursive delete @[ $$(basename $$(pwd)) = .testing ] rm -rf work results diff --git a/.travis.yml b/.travis.yml index ac7cab1b14..6bf509ce8c 100644 --- a/.travis.yml +++ b/.travis.yml @@ -5,19 +5,16 @@ language: c dist: bionic -# --depth flag is breaking our merge, try disabling it -# NOTE: We may be able to go back to depth=50 in production -git: - depth: false - addons: apt: sources: - ubuntu-toolchain-r-test packages: - - tcsh pkg-config netcdf-bin libnetcdf-dev libnetcdff-dev openmpi-bin libopenmpi-dev gfortran + - tcsh pkg-config netcdf-bin libnetcdf-dev libnetcdff-dev gfortran + - mpich libmpich-dev - doxygen graphviz flex bison cmake - python-numpy python-netcdf4 + - bc jobs: include: @@ -31,7 +28,7 @@ jobs: - test ! -s doxy_errors - env: - - JOB="Configuration testing" + - JOB="x86 Configuration testing" - DO_REGRESSION_TESTS=false - MKMF_TEMPLATE=linux-ubuntu-xenial-gnu.mk script: @@ -39,14 +36,13 @@ jobs: - echo 'Build executables...' && echo -en 'travis_fold:start:script.1\\r' - make all - echo -en 'travis_fold:end:script.1\\r' - - echo 'Running tests...' && echo -en 'travis_fold:start:script.2\\r' - - make test - - echo -en 'travis_fold:end:script.2\\r' + - make -k -s test + - make test.summary # NOTE: Code coverage upload is here to reduce load imbalance - if: type = pull_request env: - - JOB="Regression testing" + - JOB="x86 Regression testing" - DO_REGRESSION_TESTS=true - REPORT_COVERAGE=true - MKMF_TEMPLATE=linux-ubuntu-xenial-gnu.mk @@ -57,6 +53,19 @@ jobs: - echo 'Build executables...' && echo -en 'travis_fold:start:script.1\\r' - make build.regressions - echo -en 'travis_fold:end:script.1\\r' - - echo 'Running tests...' && echo -en 'travis_fold:start:script.2\\r' - - make test.regressions - - echo -en 'travis_fold:end:script.2\\r' + - make -k -s test.regressions + - make test.summary + + - arch: arm64 + env: + - JOB="ARM64 Configuration testing" + - DO_REGRESSION_TESTS=false + - DO_REPRO_TESTS=false + - MKMF_TEMPLATE=linux-ubuntu-xenial-gnu.mk + script: + - cd .testing + - echo 'Build executables...' && echo -en 'travis_fold:start:script.1\\r' + - make all + - echo -en 'travis_fold:end:script.1\\r' + - make -k -s test + - make test.summary diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 9bc71dd15f..ed6e66e0ae 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -1353,7 +1353,7 @@ subroutine build_rho_grid( G, GV, US, h, tv, dzInterface, remapCS, CS ) integer :: nz integer :: i, j, k real :: nominalDepth ! Depth of the bottom of the ocean, positive downward [H ~> m or kg m-2] - real, dimension(SZK_(GV)+1) :: zOld, zNew + real, dimension(SZK_(GV)+1) :: zOld, zNew ! Old and new interface heights [H ~> m or kg m-2] real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2] #ifdef __DO_SAFETY_CHECKS__ real :: totalThickness diff --git a/src/ALE/coord_adapt.F90 b/src/ALE/coord_adapt.F90 index 42ae0ee245..8fa4b09fc5 100644 --- a/src/ALE/coord_adapt.F90 +++ b/src/ALE/coord_adapt.F90 @@ -209,7 +209,7 @@ subroutine build_adapt_column(CS, G, GV, US, tv, i, j, zInt, tInt, sInt, h, zNex ! a positive curvature means we're too light relative to adjacent columns, ! so del2sigma needs to be positive too (push the interface deeper) call calculate_density_derivs(tInt(i,j,:), sInt(i,j,:), zInt(i,j,:) * (GV%H_to_RZ * GV%g_Earth), & - alpha, beta, tv%eqn_of_state, (/1,nz/) ) !### This should be (/1,nz+1/) - see 25 lines below. + alpha, beta, tv%eqn_of_state, (/1,nz+1/) ) do K = 2, nz ! TODO make lower bound here configurable dh_d2s(K) = del2sigma(K) * (0.5 * (h(i,j,k-1) + h(i,j,k))) / & diff --git a/src/ALE/coord_rho.F90 b/src/ALE/coord_rho.F90 index 5299c74b1b..dce802ff3c 100644 --- a/src/ALE/coord_rho.F90 +++ b/src/ALE/coord_rho.F90 @@ -91,7 +91,7 @@ subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, & h_neglect, h_neglect_edge) type(rho_CS), intent(in) :: CS !< coord_rho control structure integer, intent(in) :: nz !< Number of levels on source grid (i.e. length of h, T, S) - real, intent(in) :: depth !< Depth of ocean bottom (positive in m) + real, intent(in) :: depth !< Depth of ocean bottom (positive downward) [H ~> m or kg m-2] real, dimension(nz), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] real, dimension(nz), intent(in) :: T !< Temperature for source column [degC] real, dimension(nz), intent(in) :: S !< Salinity for source column [ppt] @@ -111,7 +111,7 @@ subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, & real, dimension(nz) :: densities ! Layer density [R ~> kg m-3] real, dimension(nz+1) :: xTmp ! Temporary positions [H ~> m or kg m-2] real, dimension(CS%nk) :: h_new ! New thicknesses [H ~> m or kg m-2] - real, dimension(CS%nk+1) :: x1 + real, dimension(CS%nk+1) :: x1 ! Interface heights [H ~> m or kg m-2] ! Construct source column with vanished layers removed (stored in h_nv) call copy_finite_thicknesses(nz, h, CS%min_thickness, count_nonzero_layers, h_nv, mapping) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index aff4860a21..a044f95893 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -45,7 +45,7 @@ module MOM use MOM_spatial_means, only : global_mass_integral use MOM_time_manager, only : time_type, real_to_time, time_type_to_real, operator(+) use MOM_time_manager, only : operator(-), operator(>), operator(*), operator(/) -use MOM_time_manager, only : operator(>=), increment_date +use MOM_time_manager, only : operator(>=), operator(==), increment_date use MOM_unit_tests, only : unit_tests use coupler_types_mod, only : coupler_type_send_data, coupler_1d_bc_type, coupler_type_spawn @@ -1013,7 +1013,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & Time_local + real_to_time(US%T_to_s*(bbl_time_int-dt)), CS%diag) ! Calculate the BBL properties and store them inside visc (u,h). call cpu_clock_begin(id_clock_BBL_visc) - call set_viscous_BBL(CS%u(:,:,:), CS%v(:,:,:), CS%h, CS%tv, CS%visc, G, GV, US, & + call set_viscous_BBL(CS%u, CS%v, CS%h, CS%tv, CS%visc, G, GV, US, & CS%set_visc_CSp, symmetrize=.true.) call cpu_clock_end(id_clock_BBL_visc) if (showCallTree) call callTree_wayPoint("done with set_viscous_BBL (step_MOM)") @@ -1404,7 +1404,8 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS logical :: skip_diffusion integer :: id_eta_diff_end - integer, pointer :: accumulated_time => NULL() + type(time_type), pointer :: accumulated_time => NULL() + type(time_type), pointer :: vertical_time => NULL() integer :: i,j,k integer :: is, ie, js, je, isd, ied, jsd, jed @@ -1426,32 +1427,30 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS call cpu_clock_begin(id_clock_offline_tracer) call extract_offline_main(CS%offline_CSp, uhtr, vhtr, eatr, ebtr, h_end, accumulated_time, & - dt_offline, dt_offline_vertical, skip_diffusion) + vertical_time, dt_offline, dt_offline_vertical, skip_diffusion) Time_end = increment_date(Time_start, seconds=floor(time_interval+0.001)) call enable_averaging(time_interval, Time_end, CS%diag) ! Check to see if this is the first iteration of the offline interval - if (accumulated_time==0) then + if (accumulated_time == real_to_time(0.0)) then first_iter = .true. else ! This is probably unnecessary but is used to guard against unwanted behavior first_iter = .false. endif - ! Check to see if vertical tracer functions should be done - if ( mod(accumulated_time, floor(US%T_to_s*dt_offline_vertical + 1e-6)) == 0 ) then + ! Check to see if vertical tracer functions should be done + if (first_iter .or. (accumulated_time >= vertical_time)) then do_vertical = .true. + vertical_time = accumulated_time + real_to_time(US%T_to_s*dt_offline_vertical) else do_vertical = .false. endif ! Increment the amount of time elapsed since last read and check if it's time to roll around - accumulated_time = mod(accumulated_time + int(time_interval), floor(US%T_to_s*dt_offline+1e-6)) - if (accumulated_time==0) then - last_iter = .true. - else - last_iter = .false. - endif + accumulated_time = accumulated_time + real_to_time(time_interval) + + last_iter = (accumulated_time >= real_to_time(US%T_to_s*dt_offline)) if (CS%use_ALE_algorithm) then ! If this is the first iteration in the offline timestep, then we need to read in fields and @@ -1566,6 +1565,10 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS fluxes%fluxes_used = .true. + if (last_iter) then + accumulated_time = real_to_time(0.0) + endif + call cpu_clock_end(id_clock_offline_tracer) end subroutine step_offline @@ -2204,7 +2207,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & allocate(CS%tv%frazil(isd:ied,jsd:jed)) ; CS%tv%frazil(:,:) = 0.0 endif if (bound_salinity) then - allocate(CS%tv%salt_deficit(isd:ied,jsd:jed)) ; CS%tv%salt_deficit(:,:)=0.0 + allocate(CS%tv%salt_deficit(isd:ied,jsd:jed)) ; CS%tv%salt_deficit(:,:) = 0.0 endif if (bulkmixedlayer .or. use_temperature) then @@ -2369,20 +2372,17 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & if (associated(sponge_in_CSp)) then ! TODO: Implementation and testing of non-ALE spong rotation - call MOM_error(FATAL, "Index rotation of non-ALE sponge is not yet " & - // "implemented.") + call MOM_error(FATAL, "Index rotation of non-ALE sponge is not yet implemented.") endif if (associated(ALE_sponge_in_CSp)) then - call rotate_ALE_sponge(ALE_sponge_in_CSp, G_in, CS%ALE_sponge_CSp, G, & - turns, param_file) - call update_ALE_sponge_field(CS%ALE_sponge_CSp, T_in, CS%T) - call update_ALE_sponge_field(CS%ALE_sponge_CSp, S_in, CS%S) + call rotate_ALE_sponge(ALE_sponge_in_CSp, G_in, CS%ALE_sponge_CSp, G, turns, param_file) + call update_ALE_sponge_field(CS%ALE_sponge_CSp, T_in, G, GV, CS%T) + call update_ALE_sponge_field(CS%ALE_sponge_CSp, S_in, G, GV, CS%S) endif if (associated(OBC_in)) & - call rotate_OBC_init(OBC_in, G, GV, US, param_file, CS%tv, restart_CSp, & - CS%OBC) + call rotate_OBC_init(OBC_in, G, GV, US, param_file, CS%tv, restart_CSp, CS%OBC) deallocate(u_in) deallocate(v_in) @@ -2427,7 +2427,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & G => CS%G if (CS%debug .or. CS%G%symmetric) then call clone_MOM_domain(CS%G%Domain, CS%G%Domain_aux, symmetric=.false.) - else ; CS%G%Domain_aux => CS%G%Domain ;endif + else ; CS%G%Domain_aux => CS%G%Domain ; endif G%ke = GV%ke endif diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 957290a04f..e6f11987a2 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -302,6 +302,7 @@ module MOM_barotropic integer :: id_visc_rem_u = -1, id_visc_rem_v = -1, id_eta_cor = -1 integer :: id_ubt = -1, id_vbt = -1, id_eta_bt = -1, id_ubtav = -1, id_vbtav = -1 integer :: id_ubt_st = -1, id_vbt_st = -1, id_eta_st = -1 + integer :: id_ubtdt = -1, id_vbtdt = -1 integer :: id_ubt_hifreq = -1, id_vbt_hifreq = -1, id_eta_hifreq = -1 integer :: id_uhbt_hifreq = -1, id_vhbt_hifreq = -1, id_eta_pred_hifreq = -1 integer :: id_gtotn = -1, id_gtots = -1, id_gtote = -1, id_gtotw = -1 @@ -476,10 +477,14 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, ! sums less than one due to viscous losses. Nondimensional. real, dimension(SZIB_(G),SZJ_(G)) :: & av_rem_u, & ! The weighted average of visc_rem_u, nondimensional. - tmp_u ! A temporary array at u points. + tmp_u, & ! A temporary array at u points. + ubt_st, & ! The zonal barotropic velocity at the start of timestep [L T-1 ~> m s-1]. + ubt_dt ! The zonal barotropic velocity tendency [L T-2 ~> m s-2]. real, dimension(SZI_(G),SZJB_(G)) :: & av_rem_v, & ! The weighted average of visc_rem_v, nondimensional. - tmp_v ! A temporary array at v points. + tmp_v, & ! A temporary array at v points. + vbt_st, & ! The meridional barotropic velocity at the start of timestep [L T-1 ~> m s-1]. + vbt_dt ! The meridional barotropic velocity tendency [L T-2 ~> m s-2]. real, dimension(SZI_(G),SZJ_(G)) :: & e_anom ! The anomaly in the sea surface height or column mass ! averaged between the beginning and end of the time step, @@ -1547,6 +1552,17 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, haloshift=1, scalar_pair=.true.) endif + if (CS%id_ubtdt > 0) then + do j=js-1,je+1 ; do I=is-1,ie + ubt_st(I,j) = ubt(I,j) + enddo ; enddo + endif + if (CS%id_vbtdt > 0) then + do J=js-1,je ; do i=is-1,ie+1 + vbt_st(i,J) = vbt(i,J) + enddo ; enddo + endif + if (query_averaging_enabled(CS%diag)) then if (CS%id_eta_st > 0) call post_data(CS%id_eta_st, eta(isd:ied,jsd:jed), CS%diag) if (CS%id_ubt_st > 0) call post_data(CS%id_ubt_st, ubt(IsdB:IedB,jsd:jed), CS%diag) @@ -2300,6 +2316,19 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, enddo ; enddo call post_data(CS%id_Corv_bt, Corv_bt_sum(isd:ied,JsdB:JedB), CS%diag) endif + if (CS%id_ubtdt > 0) then + do j=js,je ; do I=is-1,ie + ubt_dt(I,j) = (ubt_wtd(I,j) - ubt_st(I,j))*Idt + enddo ; enddo + call post_data(CS%id_ubtdt, ubt_dt(IsdB:IedB,jsd:jed), CS%diag) + endif + if (CS%id_vbtdt > 0) then + do J=js-1,je ; do i=is,ie + vbt_dt(i,J) = (vbt_wtd(i,J) - vbt_st(i,J))*Idt + enddo ; enddo + call post_data(CS%id_vbtdt, vbt_dt(isd:ied,JsdB:JedB), CS%diag) + endif + if (CS%id_ubtforce > 0) call post_data(CS%id_ubtforce, BT_force_u(IsdB:IedB,jsd:jed), CS%diag) if (CS%id_vbtforce > 0) call post_data(CS%id_vbtforce, BT_force_v(isd:ied,JsdB:JedB), CS%diag) if (CS%id_uaccel > 0) call post_data(CS%id_uaccel, u_accel_bt(IsdB:IedB,jsd:jed), CS%diag) @@ -4379,6 +4408,10 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, 'Barotropic zonal acceleration from baroclinic terms', 'm s-2', conversion=US%L_T2_to_m_s2) CS%id_vbtforce = register_diag_field('ocean_model', 'vbtforce', diag%axesCv1, Time, & 'Barotropic meridional acceleration from baroclinic terms', 'm s-2', conversion=US%L_T2_to_m_s2) + CS%id_ubtdt = register_diag_field('ocean_model', 'ubt_dt', diag%axesCu1, Time, & + 'Barotropic zonal acceleration', 'm s-2', conversion=US%L_T2_to_m_s2) + CS%id_vbtdt = register_diag_field('ocean_model', 'vbt_dt', diag%axesCv1, Time, & + 'Barotropic meridional acceleration', 'm s-2', conversion=US%L_T2_to_m_s2) CS%id_eta_bt = register_diag_field('ocean_model', 'eta_bt', diag%axesT1, Time, & 'Barotropic end SSH', thickness_units, conversion=GV%H_to_m) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 2517846b4c..5a20e60b04 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -388,10 +388,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s !--- begin set up for group halo pass cont_stencil = continuity_stencil(CS%continuity_CSp) - !### Apart from circle_OBCs halo for eta could be 1, but halo>=3 is required - !### to match circle_OBCs solutions. Why? call cpu_clock_begin(id_clock_pass) - call create_group_pass(CS%pass_eta, eta, G%Domain) !### , halo=1) + call create_group_pass(CS%pass_eta, eta, G%Domain, halo=1) call create_group_pass(CS%pass_visc_rem, CS%visc_rem_u, CS%visc_rem_v, G%Domain, & To_All+SCALAR_PAIR, CGRID_NE, halo=max(1,cont_stencil)) call create_group_pass(CS%pass_uvp, up, vp, G%Domain, halo=max(1,cont_stencil)) diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index 8c6e7d4299..d6a5186be3 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -120,6 +120,10 @@ module MOM_dynamics_unsplit real, pointer, dimension(:,:) :: tauy_bot => NULL() !< frictional y-bottom stress from the ocean !! to the seafloor [R L Z T-2 ~> Pa] + logical :: use_correct_dt_visc !< If true, use the correct timestep in the viscous terms applied + !! in the first predictor step with the unsplit time stepping scheme, + !! and in the calculation of the turbulent mixed layer properties + !! for viscosity. The default should be true, but it is false. logical :: debug !< If true, write verbose checksums for debugging purposes. logical :: module_is_initialized = .false. !< Record whether this mouled has been initialzed. @@ -228,6 +232,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vp, vpp ! Predicted meridional velocities [L T-1 ~> m s-1] real, dimension(:,:), pointer :: p_surf => NULL() real :: dt_pred ! The time step for the predictor part of the baroclinic time stepping [T ~> s]. + real :: dt_visc ! The time step for a part of the update due to viscosity [T ~> s]. logical :: dyn_p_surf integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke @@ -255,8 +260,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & ! diffu = horizontal viscosity terms (u,h) call enable_averages(dt, Time_local, CS%diag) call cpu_clock_begin(id_clock_horvisc) - call horizontal_viscosity(u, v, h, CS%diffu, CS%diffv, MEKE, Varmix, & - G, GV, US, CS%hor_visc_CSp) + call horizontal_viscosity(u, v, h, CS%diffu, CS%diffv, MEKE, Varmix, G, GV, US, CS%hor_visc_CSp) call cpu_clock_end(id_clock_horvisc) call disable_averaging(CS%diag) @@ -323,31 +327,29 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & ! up = u + dt_pred * (PFu + CAu) call cpu_clock_begin(id_clock_mom_update) do k=1,nz ; do j=js,je ; do I=Isq,Ieq - up(I,j,k) = G%mask2dCu(I,j) * (u(I,j,k) + dt_pred * & - (CS%PFu(I,j,k) + CS%CAu(I,j,k))) + up(I,j,k) = G%mask2dCu(I,j) * (u(I,j,k) + dt_pred * (CS%PFu(I,j,k) + CS%CAu(I,j,k))) enddo ; enddo ; enddo do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie - vp(i,J,k) = G%mask2dCv(i,J) * (v(i,J,k) + dt_pred * & - (CS%PFv(i,J,k) + CS%CAv(i,J,k))) + vp(i,J,k) = G%mask2dCv(i,J) * (v(i,J,k) + dt_pred * (CS%PFv(i,J,k) + CS%CAv(i,J,k))) enddo ; enddo ; enddo call cpu_clock_end(id_clock_mom_update) if (CS%debug) then call MOM_state_chksum("Predictor 1", up, vp, h_av, uh, vh, G, GV, US) - call MOM_accel_chksum("Predictor 1 accel", CS%CAu, CS%CAv, CS%PFu, CS%PFv,& + call MOM_accel_chksum("Predictor 1 accel", CS%CAu, CS%CAv, CS%PFu, CS%PFv, & CS%diffu, CS%diffv, G, GV, US) endif ! up <- up + dt/2 d/dz visc d/dz up call cpu_clock_begin(id_clock_vertvisc) call enable_averages(dt, Time_local, CS%diag) - call set_viscous_ML(u, v, h_av, tv, forces, visc, dt*0.5, G, GV, US, & - CS%set_visc_CSp) + dt_visc = 0.5*dt ; if (CS%use_correct_dt_visc) dt_visc = dt + call set_viscous_ML(u, v, h_av, tv, forces, visc, dt_visc, G, GV, US, CS%set_visc_CSp) call disable_averaging(CS%diag) - !### I think that the time steps in the next two calls should be dt_pred. - call vertvisc_coef(up, vp, h_av, forces, visc, dt*0.5, G, GV, US, & - CS%vertvisc_CSp, CS%OBC) - call vertvisc(up, vp, h_av, forces, visc, dt*0.5, CS%OBC, CS%ADp, CS%CDp, & + + dt_visc = 0.5*dt ; if (CS%use_correct_dt_visc) dt_visc = dt_pred + call vertvisc_coef(up, vp, h_av, forces, visc, dt_visc, G, GV, US, CS%vertvisc_CSp, CS%OBC) + call vertvisc(up, vp, h_av, forces, visc, dt_visc, CS%OBC, CS%ADp, CS%CDp, & G, GV, US, CS%vertvisc_CSp, Waves=Waves) call cpu_clock_end(id_clock_vertvisc) call pass_vector(up, vp, G%Domain, clock=id_clock_pass) @@ -355,8 +357,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & ! uh = up * hp ! h_av = hp + dt/2 div . uh call cpu_clock_begin(id_clock_continuity) - call continuity(up, vp, hp, h_av, uh, vh, (0.5*dt), G, GV, US, & - CS%continuity_CSp, OBC=CS%OBC) + call continuity(up, vp, hp, h_av, uh, vh, (0.5*dt), G, GV, US, CS%continuity_CSp, OBC=CS%OBC) call cpu_clock_end(id_clock_continuity) call pass_var(h_av, G%Domain, clock=id_clock_pass) call pass_vector(uh, vh, G%Domain, clock=id_clock_pass) @@ -392,25 +393,22 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & ! upp = u + dt/2 * ( PFu + CAu ) call cpu_clock_begin(id_clock_mom_update) do k=1,nz ; do j=js,je ; do I=Isq,Ieq - upp(I,j,k) = G%mask2dCu(I,j) * (u(I,j,k) + dt * 0.5 * & - (CS%PFu(I,j,k) + CS%CAu(I,j,k))) + upp(I,j,k) = G%mask2dCu(I,j) * (u(I,j,k) + dt * 0.5 * (CS%PFu(I,j,k) + CS%CAu(I,j,k))) enddo ; enddo ; enddo do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie - vpp(i,J,k) = G%mask2dCv(i,J) * (v(i,J,k) + dt * 0.5 * & - (CS%PFv(i,J,k) + CS%CAv(i,J,k))) + vpp(i,J,k) = G%mask2dCv(i,J) * (v(i,J,k) + dt * 0.5 * (CS%PFv(i,J,k) + CS%CAv(i,J,k))) enddo ; enddo ; enddo call cpu_clock_end(id_clock_mom_update) if (CS%debug) then call MOM_state_chksum("Predictor 2", upp, vpp, h_av, uh, vh, G, GV, US) - call MOM_accel_chksum("Predictor 2 accel", CS%CAu, CS%CAv, CS%PFu, CS%PFv,& + call MOM_accel_chksum("Predictor 2 accel", CS%CAu, CS%CAv, CS%PFu, CS%PFv, & CS%diffu, CS%diffv, G, GV, US) endif ! upp <- upp + dt/2 d/dz visc d/dz upp call cpu_clock_begin(id_clock_vertvisc) - call vertvisc_coef(upp, vpp, hp, forces, visc, dt*0.5, G, GV, US, & - CS%vertvisc_CSp, CS%OBC) + call vertvisc_coef(upp, vpp, hp, forces, visc, dt*0.5, G, GV, US, CS%vertvisc_CSp, CS%OBC) call vertvisc(upp, vpp, hp, forces, visc, dt*0.5, CS%OBC, CS%ADp, CS%CDp, & G, GV, US, CS%vertvisc_CSp, Waves=Waves) call cpu_clock_end(id_clock_vertvisc) @@ -419,8 +417,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & ! uh = upp * hp ! h = hp + dt/2 div . uh call cpu_clock_begin(id_clock_continuity) - call continuity(upp, vpp, hp, h, uh, vh, (dt*0.5), G, GV, US, & - CS%continuity_CSp, OBC=CS%OBC) + call continuity(upp, vpp, hp, h, uh, vh, (dt*0.5), G, GV, US, CS%continuity_CSp, OBC=CS%OBC) call cpu_clock_end(id_clock_continuity) call pass_var(h, G%Domain, clock=id_clock_pass) call pass_vector(uh, vh, G%Domain, clock=id_clock_pass) @@ -470,12 +467,10 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & call open_boundary_zero_normal_flow(CS%OBC, G, CS%CAu, CS%CAv) endif do k=1,nz ; do j=js,je ; do I=Isq,Ieq - u(I,j,k) = G%mask2dCu(I,j) * (u(I,j,k) + dt * & - (CS%PFu(I,j,k) + CS%CAu(I,j,k))) + u(I,j,k) = G%mask2dCu(I,j) * (u(I,j,k) + dt * (CS%PFu(I,j,k) + CS%CAu(I,j,k))) enddo ; enddo ; enddo do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie - v(i,J,k) = G%mask2dCv(i,J) * (v(i,J,k) + dt * & - (CS%PFv(i,J,k) + CS%CAv(i,J,k))) + v(i,J,k) = G%mask2dCv(i,J) * (v(i,J,k) + dt * (CS%PFv(i,J,k) + CS%CAv(i,J,k))) enddo ; enddo ; enddo ! u <- u + dt d/dz visc d/dz u @@ -634,6 +629,11 @@ subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, US, param_file, diag, CS CS%diag => diag + call get_param(param_file, mdl, "FIX_UNSPLIT_DT_VISC_BUG", CS%use_correct_dt_visc, & + "If true, use the correct timestep in the viscous terms applied in the first "//& + "predictor step with the unsplit time stepping scheme, and in the calculation "//& + "of the turbulent mixed layer properties for viscosity with unsplit or "//& + "unsplit_RK2. The default should be true.", default=.false.) call get_param(param_file, mdl, "DEBUG", CS%debug, & "If true, write out verbose debugging data.", & default=.false., debuggingParam=.true.) diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index d3adfaa194..e3ec48ff58 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -123,6 +123,9 @@ module MOM_dynamics_unsplit_RK2 !! the extent to which the treatment of gravity waves !! is forward-backward (0) or simulated backward !! Euler (1). 0 is almost always used. + logical :: use_correct_dt_visc !< If true, use the correct timestep in the calculation of the + !! turbulent mixed layer properties for viscosity. + !! The default should be true, but it is false. logical :: debug !< If true, write verbose checksums for debugging purposes. logical :: module_is_initialized = .false. !< Record whether this mouled has been initialzed. @@ -238,8 +241,8 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: up ! Predicted zonal velocities [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vp ! Predicted meridional velocities [L T-1 ~> m s-1] real, dimension(:,:), pointer :: p_surf => NULL() - real :: dt_pred ! The time step for the predictor part of the baroclinic - ! time stepping [T ~> s]. + real :: dt_pred ! The time step for the predictor part of the baroclinic time stepping [T ~> s] + real :: dt_visc ! The time step for a part of the update due to viscosity [T ~> s] logical :: dyn_p_surf integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke @@ -280,17 +283,15 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, call cpu_clock_begin(id_clock_continuity) ! This is a duplicate calculation of the last continuity from the previous step ! and could/should be optimized out. -AJA - call continuity(u_in, v_in, h_in, hp, uh, vh, dt_pred, G, GV, US, & - CS%continuity_CSp, OBC=CS%OBC) + call continuity(u_in, v_in, h_in, hp, uh, vh, dt_pred, G, GV, US, CS%continuity_CSp, OBC=CS%OBC) call cpu_clock_end(id_clock_continuity) call pass_var(hp, G%Domain, clock=id_clock_pass) call pass_vector(uh, vh, G%Domain, clock=id_clock_pass) ! h_av = (h + hp)/2 (used in PV denominator) call cpu_clock_begin(id_clock_mom_update) - do k=1,nz - do j=js-2,je+2 ; do i=is-2,ie+2 - h_av(i,j,k) = (h_in(i,j,k) + hp(i,j,k)) * 0.5 + do k=1,nz ; do j=js-2,je+2 ; do i=is-2,ie+2 + h_av(i,j,k) = (h_in(i,j,k) + hp(i,j,k)) * 0.5 enddo ; enddo ; enddo call cpu_clock_end(id_clock_mom_update) @@ -305,8 +306,7 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, if (dyn_p_surf) then ; do j=js-2,je+2 ; do i=is-2,ie+2 p_surf(i,j) = 0.5*p_surf_begin(i,j) + 0.5*p_surf_end(i,j) enddo ; enddo ; endif - call PressureForce(h_in, tv, CS%PFu, CS%PFv, G, GV, US, & - CS%PressureForce_CSp, CS%ALE_CSp, p_surf) + call PressureForce(h_in, tv, CS%PFu, CS%PFv, G, GV, US, CS%PressureForce_CSp, CS%ALE_CSp, p_surf) call cpu_clock_end(id_clock_pres) call pass_vector(CS%PFu, CS%PFv, G%Domain, clock=id_clock_pass) call pass_vector(CS%CAu, CS%CAv, G%Domain, clock=id_clock_pass) @@ -339,11 +339,11 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, ! up[n-1/2] <- up*[n-1/2] + dt/2 d/dz visc d/dz up[n-1/2] call cpu_clock_begin(id_clock_vertvisc) call enable_averages(dt, Time_local, CS%diag) - call set_viscous_ML(up, vp, h_av, tv, forces, visc, dt_pred, G, GV, US, & - CS%set_visc_CSp) + dt_visc = dt_pred ; if (CS%use_correct_dt_visc) dt_visc = dt + call set_viscous_ML(up, vp, h_av, tv, forces, visc, dt_visc, G, GV, US, CS%set_visc_CSp) call disable_averaging(CS%diag) - call vertvisc_coef(up, vp, h_av, forces, visc, dt_pred, G, GV, US, & - CS%vertvisc_CSp, CS%OBC) + + call vertvisc_coef(up, vp, h_av, forces, visc, dt_pred, G, GV, US, CS%vertvisc_CSp, CS%OBC) call vertvisc(up, vp, h_av, forces, visc, dt_pred, CS%OBC, CS%ADp, CS%CDp, & G, GV, US, CS%vertvisc_CSp) call cpu_clock_end(id_clock_vertvisc) @@ -394,12 +394,10 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, ! up[n] <- up* + dt d/dz visc d/dz up ! u[n] <- u*[n] + dt d/dz visc d/dz u[n] call cpu_clock_begin(id_clock_vertvisc) - call vertvisc_coef(up, vp, h_av, forces, visc, dt, G, GV, US, & - CS%vertvisc_CSp, CS%OBC) + call vertvisc_coef(up, vp, h_av, forces, visc, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC) call vertvisc(up, vp, h_av, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp, & G, GV, US, CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot) - call vertvisc_coef(u_in, v_in, h_av, forces, visc, dt, G, GV, US, & - CS%vertvisc_CSp, CS%OBC) + call vertvisc_coef(u_in, v_in, h_av, forces, visc, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC) call vertvisc(u_in, v_in, h_av, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp,& G, GV, US, CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot) call cpu_clock_end(id_clock_vertvisc) @@ -593,6 +591,11 @@ subroutine initialize_dyn_unsplit_RK2(u, v, h, Time, G, GV, US, param_file, diag "If SPLIT is false and USE_RK2 is true, BEGW can be "//& "between 0 and 0.5 to damp gravity waves.", & units="nondim", default=0.0) + call get_param(param_file, mdl, "FIX_UNSPLIT_DT_VISC_BUG", CS%use_correct_dt_visc, & + "If true, use the correct timestep in the viscous terms applied in the first "//& + "predictor step with the unsplit time stepping scheme, and in the calculation "//& + "of the turbulent mixed layer properties for viscosity with unsplit or "//& + "unsplit_RK2. The default should be true.", default=.false.) call get_param(param_file, mdl, "DEBUG", CS%debug, & "If true, write out verbose debugging data.", & default=.false., debuggingParam=.true.) diff --git a/src/core/MOM_isopycnal_slopes.F90 b/src/core/MOM_isopycnal_slopes.F90 index 53f0d59294..fa60fb821d 100644 --- a/src/core/MOM_isopycnal_slopes.F90 +++ b/src/core/MOM_isopycnal_slopes.F90 @@ -242,7 +242,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & slope_x(I,j,K) = 0.0 endif - if (present_N2_u) N2_u(I,j,k) = G_Rho0 * drdz * G%mask2dCu(I,j) ! Square of Brunt-Vaisala frequency [s-2] + if (present_N2_u) N2_u(I,j,k) = G_Rho0 * drdz * G%mask2dCu(I,j) ! Square of buoyancy frequency [T-2 ~> s-2] else ! With .not.use_EOS, the layers are constant density. slope_x(I,j,K) = (Z_to_L*(e(i,j,K)-e(i+1,j,K))) * G%IdxCu(I,j) @@ -328,7 +328,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & slope_y(i,J,K) = 0.0 endif - if (present_N2_v) N2_v(i,J,k) = G_Rho0 * drdz * G%mask2dCv(i,J) ! Square of Brunt-Vaisala frequency [s-2] + if (present_N2_v) N2_v(i,J,k) = G_Rho0 * drdz * G%mask2dCv(i,J) ! Square of buoyancy frequency [T-2 ~> s-2] else ! With .not.use_EOS, the layers are constant density. slope_y(i,J,K) = (Z_to_L*(e(i,j,K)-e(i,j+1,K))) * G%IdyCv(i,J) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index ffede1c0c2..5b6dc168f4 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -265,10 +265,8 @@ module MOM_open_boundary !! velocities (or speed of characteristics) at the !! new time level (1) or the running mean (0) for velocities. !! Valid values range from 0 to 1, with a default of 0.3. - real :: rx_max !< The maximum magnitude of the baroclinic radiation - !! velocity (or speed of characteristics) [m s-1]. The - !! default value is 10 m s-1. - !### The description above seems inconsistent with the code, and the units should be [nondim]. + real :: rx_max !< The maximum magnitude of the baroclinic radiation velocity (or speed of + !! characteristics) in units of grid points per timestep [nondim]. logical :: OBC_pe !< Is there an open boundary on this tile? type(remapping_CS), pointer :: remap_CS !< ALE remapping control structure for segments only type(OBC_registry_type), pointer :: OBC_Reg => NULL() !< Registry type for boundaries @@ -500,13 +498,11 @@ subroutine open_boundary_config(G, US, param_file, OBC) call initialize_segment_data(G, OBC, param_file) if (open_boundary_query(OBC, apply_open_OBC=.true.)) then - !### I think that OBC%rx_max as used is actually nondimensional, with effective - ! units of grid points per time step. call get_param(param_file, mdl, "OBC_RADIATION_MAX", OBC%rx_max, & - "The maximum magnitude of the baroclinic radiation "//& - "velocity (or speed of characteristics). This is only "//& + "The maximum magnitude of the baroclinic radiation velocity (or speed of "//& + "characteristics), in gridpoints per timestep. This is only "//& "used if one of the open boundary segments is using Orlanski.", & - units="m s-1", default=10.0) !### Should the units here be "nondim"? + units="nondim", default=10.0) !### Should the default be changed to 1.0? call get_param(param_file, mdl, "OBC_RAD_VEL_WT", OBC%gamma_uv, & "The relative weighting for the baroclinic radiation "//& "velocities (or speed of characteristics) at the new "//& @@ -1804,9 +1800,9 @@ subroutine open_boundary_impose_land_mask(OBC, G, areaCu, areaCv, US) I=segment%HI%IsdB do j=segment%HI%jsd,segment%HI%jed if (segment%direction == OBC_DIRECTION_E) then - areaCu(I,j) = G%areaT(i,j) ! Both of these are in [L2] + areaCu(I,j) = G%areaT(i,j) ! Both of these are in [L2 ~> m2] else ! West - areaCu(I,j) = G%areaT(i+1,j) ! Both of these are in [L2] + areaCu(I,j) = G%areaT(i+1,j) ! Both of these are in [L2 ~> m2] endif enddo else @@ -1814,9 +1810,9 @@ subroutine open_boundary_impose_land_mask(OBC, G, areaCu, areaCv, US) J=segment%HI%JsdB do i=segment%HI%isd,segment%HI%ied if (segment%direction == OBC_DIRECTION_S) then - areaCv(i,J) = G%areaT(i,j+1) ! Both of these are in [L2] + areaCv(i,J) = G%areaT(i,j+1) ! Both of these are in [L2 ~> m2] else ! North - areaCu(i,J) = G%areaT(i,j) ! Both of these are in [L2] + areaCu(i,J) = G%areaT(i,j) ! Both of these are in [L2 ~> m2] endif enddo endif @@ -1906,7 +1902,7 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, US, dt) real :: dhdt, dhdx, dhdy ! One-point differences in time or space [L T-1 ~> m s-1] real :: gamma_u, gamma_2 ! Fractional weightings of new values [nondim] real :: tau ! A local nudging timescale [T ~> s] - real :: rx_max, ry_max ! coefficients for radiation [nondim] or [L2 T-2 ~> m2 s-2] + real :: rx_max, ry_max ! coefficients for radiation [nondim] real :: rx_new, rx_avg ! coefficients for radiation [nondim] or [L2 T-2 ~> m2 s-2] real :: ry_new, ry_avg ! coefficients for radiation [nondim] or [L2 T-2 ~> m2 s-2] real :: cff_new, cff_avg ! denominator in oblique [L2 T-2 ~> m2 s-2] diff --git a/src/diagnostics/MOM_wave_speed.F90 b/src/diagnostics/MOM_wave_speed.F90 index 65a23e0fa2..9da2963c16 100644 --- a/src/diagnostics/MOM_wave_speed.F90 +++ b/src/diagnostics/MOM_wave_speed.F90 @@ -1246,7 +1246,7 @@ subroutine tridiag_det(a, b, c, nrows, lam, det_out, ddet_out, row_scale) rscl = 1.0 ; if (present(row_scale)) rscl = row_scale det(1) = 1.0 ; ddet(1) = 0.0 - det(2) = b(2)-lam ; ddet(2) = -1.0 + if (nrows > 1) then ; det(2) = b(2)-lam ; ddet(2) = -1.0 ; endif do n=3,nrows det(n) = rscl*(b(n)-lam)*det(n-1) - rscl*(a(n)*c(n-1))*det(n-2) ddet(n) = rscl*(b(n)-lam)*ddet(n-1) - rscl*(a(n)*c(n-1))*ddet(n-2) - det(n-1) diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index 0fe937a173..4e12abaa5b 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -31,11 +31,11 @@ ! ! For interpolation between h and u grids, we use the following relations: ! -! h->u: f_u[ig] = 0.5 * (f_h[ ig ] + f_h[ig+1]) -! f_u[i1] = 0.5 * (f_h[i1-1] + f_h[ i1 ]) +! h->u: f_u(ig) = 0.5 * (f_h( ig ) + f_h(ig+1)) +! f_u(i1) = 0.5 * (f_h(i1-1) + f_h( i1 )) ! -! u->h: f_h[ig] = 0.5 * (f_u[ig-1] + f_u[ ig ]) -! f_h[i1] = 0.5 * (f_u[ i1 ] + f_u[i1+1]) +! u->h: f_h(ig) = 0.5 * (f_u(ig-1) + f_u( ig )) +! f_h(i1) = 0.5 * (f_u( i1 ) + f_u(i1+1)) ! ! where ig is the grid index and i1 is the 1-based index. That is, a 1-based ! u-point is ahead of its matching h-point in non-symmetric mode, but behind @@ -110,11 +110,11 @@ module MOM_diag_remap character(len=16) :: diag_coord_name = '' !< A name for the purpose of run-time parameters character(len=8) :: diag_module_suffix = '' !< The suffix for the module to appear in diag_table type(remapping_CS) :: remap_cs !< Remapping control structure use for this axes - type(regridding_CS) :: regrid_cs !< Regridding control structure that defines the coordiantes for this axes + type(regridding_CS) :: regrid_cs !< Regridding control structure that defines the coordinates for this axes integer :: nz = 0 !< Number of vertical levels used for remapping - real, dimension(:,:,:), allocatable :: h !< Remap grid thicknesses - real, dimension(:,:,:), allocatable :: h_extensive !< Remap grid thicknesses for extensive variables - real, dimension(:), allocatable :: dz !< Nominal layer thicknesses + real, dimension(:,:,:), allocatable :: h !< Remap grid thicknesses [H ~> m or kg m-2] + real, dimension(:,:,:), allocatable :: h_extensive !< Remap grid thicknesses for extensive + !! variables [H ~> m or kg m-2] integer :: interface_axes_id = 0 !< Vertical axes id for remapping at interfaces integer :: layer_axes_id = 0 !< Vertical axes id for remapping on layers logical :: answers_2018 !< If true, use the order of arithmetic and expressions for remapping @@ -151,7 +151,6 @@ subroutine diag_remap_end(remap_cs) type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diag remapping control structure if (allocated(remap_cs%h)) deallocate(remap_cs%h) - if (allocated(remap_cs%dz)) deallocate(remap_cs%dz) remap_cs%configured = .false. remap_cs%initialized = .false. remap_cs%used = .false. @@ -274,15 +273,15 @@ function diag_remap_axes_configured(remap_cs) !! coordinates then technically we should also regenerate the !! target grid whenever T/S change. subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_target) - type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diagnostic coordinate control structure - type(ocean_grid_type), pointer :: G !< The ocean's grid type - type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - real, dimension(:, :, :), intent(in) :: h !< New thickness [H ~> m or kg m-2] - real, dimension(:, :, :), intent(in) :: T !< New temperatures [degC] - real, dimension(:, :, :), intent(in) :: S !< New salinities [ppt] - type(EOS_type), pointer :: eqn_of_state !< A pointer to the equation of state - real, dimension(:, :, :), intent(inout) :: h_target !< Where to store the new diagnostic array + type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diagnostic coordinate control structure + type(ocean_grid_type), pointer :: G !< The ocean's grid type + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(:,:,:), intent(in) :: h !< New thickness [H ~> m or kg m-2] + real, dimension(:,:,:), intent(in) :: T !< New temperatures [degC] + real, dimension(:,:,:), intent(in) :: S !< New salinities [ppt] + type(EOS_type), pointer :: eqn_of_state !< A pointer to the equation of state + real, dimension(:,:,:), intent(inout) :: h_target !< The new diagnostic thicknesses [H ~> m or kg m-2] ! Local variables real, dimension(remap_cs%nz + 1) :: zInterfaces ! Interface positions [H ~> m or kg m-2] @@ -328,9 +327,8 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_targe call build_sigma_column(get_sigma_CS(remap_cs%regrid_cs), & GV%Z_to_H*G%bathyT(i,j), sum(h(i,j,:)), zInterfaces) elseif (remap_cs%vertical_coord == coordinateMode('RHO')) then -!### I think that the conversion factor in the 2nd line should be GV%Z_to_H call build_rho_column(get_rho_CS(remap_cs%regrid_cs), G%ke, & - US%Z_to_m*G%bathyT(i,j), h(i,j,:), T(i,j,:), S(i,j,:), & + GV%Z_to_H*G%bathyT(i,j), h(i,j,:), T(i,j,:), S(i,j,:), & eqn_of_state, zInterfaces, h_neglect, h_neglect_edge) elseif (remap_cs%vertical_coord == coordinateMode('SLIGHT')) then ! call build_slight_column(remap_cs%regrid_cs,remap_cs%remap_cs, nz, & @@ -354,16 +352,16 @@ subroutine diag_remap_do_remap(remap_cs, G, GV, h, staggered_in_x, staggered_in_ type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coodinate control structure type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure - real, dimension(:,:,:), intent(in) :: h !< The current thicknesses + real, dimension(:,:,:), intent(in) :: h !< The current thicknesses [H ~> m or kg m-2] logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points - real, dimension(:,:,:), pointer :: mask !< A mask for the field - real, dimension(:,:,:), intent(in) :: field(:,:,:) !< The diagnostic field to be remapped - real, dimension(:,:,:), intent(inout) :: remapped_field !< Field remapped to new coordinate + real, dimension(:,:,:), pointer :: mask !< A mask for the field [nondim] + real, dimension(:,:,:), intent(in) :: field(:,:,:) !< The diagnostic field to be remapped [A] + real, dimension(:,:,:), intent(inout) :: remapped_field !< Field remapped to new coordinate [A] ! Local variables - real, dimension(remap_cs%nz) :: h_dest - real, dimension(size(h,3)) :: h_src - real :: h_neglect, h_neglect_edge + real, dimension(remap_cs%nz) :: h_dest ! Destination thicknesses [H ~> m or kg m-2] + real, dimension(size(h,3)) :: h_src ! A column of source thicknesses [H ~> m or kg m-2] + real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2] integer :: nz_src, nz_dest integer :: i, j, k !< Grid index integer :: i1, j1 !< 1-based index @@ -446,14 +444,15 @@ end subroutine diag_remap_do_remap !> Calculate masks for target grid subroutine diag_remap_calc_hmask(remap_cs, G, mask) - type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coodinate control structure - type(ocean_grid_type), intent(in) :: G !< Ocean grid structure - real, dimension(:,:,:), intent(out) :: mask !< h-point mask for target grid + type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coodinate control structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + real, dimension(:,:,:), intent(out) :: mask !< h-point mask for target grid [nondim] ! Local variables - real, dimension(remap_cs%nz) :: h_dest + real, dimension(remap_cs%nz) :: h_dest ! Destination thicknesses [H ~> m or kg m-2] integer :: i, j, k logical :: mask_vanished_layers - real :: h_tot, h_err + real :: h_tot ! Sum of all thicknesses [H ~> m or kg m-2] + real :: h_err ! An estimate of a negligible thickness [H ~> m or kg m-2] call assert(remap_cs%initialized, 'diag_remap_calc_hmask: remap_cs not initialized.') @@ -492,16 +491,16 @@ subroutine vertically_reintegrate_diag_field(remap_cs, G, h, h_target, staggered mask, field, reintegrated_field) type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coodinate control structure type(ocean_grid_type), intent(in) :: G !< Ocean grid structure - real, dimension(:,:,:), intent(in) :: h !< The thicknesses of the source grid - real, dimension(:,:,:), intent(in) :: h_target !< The thicknesses of the target grid + real, dimension(:,:,:), intent(in) :: h !< The thicknesses of the source grid [H ~> m or kg m-2] + real, dimension(:,:,:), intent(in) :: h_target !< The thicknesses of the target grid [H ~> m or kg m-2] logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points - real, dimension(:,:,:), pointer :: mask !< A mask for the field - real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped - real, dimension(:,:,:), intent(inout) :: reintegrated_field !< Field argument remapped to alternative coordinate + real, dimension(:,:,:), pointer :: mask !< A mask for the field [nondim] + real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped [A] + real, dimension(:,:,:), intent(inout) :: reintegrated_field !< Field argument remapped to alternative coordinate [A] ! Local variables - real, dimension(remap_cs%nz) :: h_dest - real, dimension(size(h,3)) :: h_src + real, dimension(remap_cs%nz) :: h_dest ! Destination thicknesses [H ~> m or kg m-2] + real, dimension(size(h,3)) :: h_src ! A column of source thicknesses [H ~> m or kg m-2] integer :: nz_src, nz_dest integer :: i, j, k !< Grid index integer :: i1, j1 !< 1-based index @@ -572,16 +571,16 @@ end subroutine vertically_reintegrate_diag_field subroutine vertically_interpolate_diag_field(remap_cs, G, h, staggered_in_x, staggered_in_y, & mask, field, interpolated_field) type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coodinate control structure - type(ocean_grid_type), intent(in) :: G !< Ocean grid structure - real, dimension(:,:,:), intent(in) :: h !< The current thicknesses + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + real, dimension(:,:,:), intent(in) :: h !< The current thicknesses [H ~> m or kg m-2] logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points - real, dimension(:,:,:), pointer :: mask !< A mask for the field - real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped - real, dimension(:,:,:), intent(inout) :: interpolated_field !< Field argument remapped to alternative coordinate + real, dimension(:,:,:), pointer :: mask !< A mask for the field [nondim] + real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped [A] + real, dimension(:,:,:), intent(inout) :: interpolated_field !< Field argument remapped to alternative coordinate [A] ! Local variables - real, dimension(remap_cs%nz) :: h_dest - real, dimension(size(h,3)) :: h_src + real, dimension(remap_cs%nz) :: h_dest ! Destination thicknesses [H ~> m or kg m-2] + real, dimension(size(h,3)) :: h_src ! A column of source thicknesses [H ~> m or kg m-2] integer :: nz_src, nz_dest integer :: i, j, k !< Grid index integer :: i1, j1 !< 1-based index @@ -656,20 +655,20 @@ subroutine horizontally_average_diag_field(G, GV, h, staggered_in_x, staggered_i averaged_mask) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean vertical grid structure - real, dimension(:,:,:), intent(in) :: h !< The current thicknesses + real, dimension(:,:,:), intent(in) :: h !< The current thicknesses [H ~> m or kg m-2] logical, intent(in) :: staggered_in_x !< True if the x-axis location is at u or q points logical, intent(in) :: staggered_in_y !< True if the y-axis location is at v or q points logical, intent(in) :: is_layer !< True if the z-axis location is at h points logical, intent(in) :: is_extensive !< True if the z-direction is spatially integrated (over layers) - real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped - real, dimension(:), intent(inout) :: averaged_field !< Field argument horizontally averaged - logical, dimension(:), intent(inout) :: averaged_mask !< Mask for horizontally averaged field + real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped [A] + real, dimension(:), intent(inout) :: averaged_field !< Field argument horizontally averaged [A] + logical, dimension(:), intent(inout) :: averaged_mask !< Mask for horizontally averaged field [nondim] ! Local variables real, dimension(G%isc:G%iec, G%jsc:G%jec, size(field,3)) :: volume, stuff real, dimension(size(field, 3)) :: vol_sum, stuff_sum ! nz+1 is needed for interface averages type(EFP_type), dimension(2*size(field,3)) :: sums_EFP ! Sums of volume or stuff by layer - real :: height + real :: height ! An average thickness attributed to an velocity point [H ~> m or kg m-2] integer :: i, j, k, nz integer :: i1, j1 !< 1-based index diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 891d6b3ea7..6b68cb3deb 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -607,15 +607,16 @@ subroutine shelf_calc_flux(sfc_state, fluxes, Time, time_step, CS, forces) ISS%tflux_ocn(i,j) = 0.0 endif -! haline_driving(:,:) = sfc_state%sss(i,j) - Sbdry(i,j) +! haline_driving(i,j) = sfc_state%sss(i,j) - Sbdry(i,j) enddo ! i-loop enddo ! j-loop - ! ISS%water_flux = net liquid water into the ocean [R Z T-1 ~> kg m-2 s-1] - fluxes%iceshelf_melt(:,:) = ISS%water_flux(:,:) * CS%flux_factor do j=js,je ; do i=is,ie + ! ISS%water_flux = net liquid water into the ocean [R Z T-1 ~> kg m-2 s-1] + fluxes%iceshelf_melt(i,j) = ISS%water_flux(i,j) * CS%flux_factor + if ((sfc_state%ocean_mass(i,j) > CS%col_mass_melt_threshold) .and. & (ISS%area_shelf_h(i,j) > 0.0) .and. (CS%isthermo)) then @@ -653,11 +654,10 @@ subroutine shelf_calc_flux(sfc_state, fluxes, Time, time_step, CS, forces) ISS%water_flux(i,j) = 0.0 fluxes%iceshelf_melt(i,j) = 0.0 endif ! area_shelf_h - enddo ; enddo ! i- and j-loops - ! mass flux [R Z L2 T-1 ~> kg s-1], part of ISOMIP diags. - mass_flux(:,:) = 0.0 - mass_flux(:,:) = ISS%water_flux(:,:) * ISS%area_shelf_h(:,:) + ! mass flux [R Z L2 T-1 ~> kg s-1], part of ISOMIP diags. + mass_flux(i,j) = ISS%water_flux(i,j) * ISS%area_shelf_h(i,j) + enddo ; enddo ! i- and j-loops if (CS%active_shelf_dynamics .or. CS%override_shelf_movement) then call cpu_clock_begin(id_clock_pass) @@ -690,7 +690,7 @@ subroutine shelf_calc_flux(sfc_state, fluxes, Time, time_step, CS, forces) ! advect the ice shelf, and advance the front. Calving will be in here somewhere as well.. ! when we decide on how to do it call update_ice_shelf(CS%dCS, ISS, G, US, US%s_to_T*time_step, Time, & - sfc_state%ocean_mass(:,:), coupled_GL) + sfc_state%ocean_mass, coupled_GL) endif @@ -1051,10 +1051,10 @@ subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes) enddo ; enddo balancing_area = global_area_integral(bal_frac, G) - if (balancing_area > 0.0) then !### Examine whether the rescaling should be inside the parenthesis. - balancing_flux = US%kg_m2s_to_RZ_T*(global_area_integral(ISS%water_flux, G, scale=US%RZ_T_to_kg_m2s, & - area=ISS%area_shelf_h) + & - delta_mass_shelf ) / balancing_area + if (balancing_area > 0.0) then + balancing_flux = ( US%kg_m2s_to_RZ_T*global_area_integral(ISS%water_flux, G, scale=US%RZ_T_to_kg_m2s, & + area=ISS%area_shelf_h) + & + delta_mass_shelf ) / balancing_area else balancing_flux = 0.0 endif @@ -1166,7 +1166,6 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl write(0,*) 'IG: ', G%isd, G%isc, G%iec, G%ied, G%jsd, G%jsc, G%jsd, G%jed endif - CS%Time = Time ! ### This might not be in the right place? CS%diag => diag ! Are we being called from the solo ice-sheet driver? When called by the ocean @@ -1735,7 +1734,9 @@ subroutine update_shelf_mass(G, US, CS, ISS, Time) call time_interp_external(CS%id_read_mass, Time, ISS%mass_shelf) ! This should only be done if time_interp_external did an update. - ISS%mass_shelf(:,:) = US%kg_m3_to_R*US%m_to_Z * ISS%mass_shelf(:,:) ! Rescale after time_interp + do j=js,je ; do i=is,ie + ISS%mass_shelf(i,j) = US%kg_m3_to_R*US%m_to_Z * ISS%mass_shelf(i,j) ! Rescale after time_interp + enddo ; enddo do j=js,je ; do i=is,ie ISS%area_shelf_h(i,j) = 0.0 diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index beeaf6e46a..07d928d76b 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -1108,11 +1108,11 @@ subroutine trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read_params) just_read = .false. ; if (present(just_read_params)) just_read = just_read_params call get_param(PF, mdl, "SURFACE_PRESSURE_FILE", p_surf_file, & - "The initial condition file for the surface height.", & + "The initial condition file for the surface pressure exerted by ice.", & fail_if_missing=.not.just_read, do_not_log=just_read) call get_param(PF, mdl, "SURFACE_PRESSURE_VAR", p_surf_var, & - "The initial condition variable for the surface height.", & - units="kg m-2", default="", do_not_log=just_read) !### The units here should be Pa? + "The initial condition variable for the surface pressure exerted by ice.", & + units="Pa", default="", do_not_log=just_read) call get_param(PF, mdl, "INPUTDIR", inputdir, default=".", do_not_log=.true.) filename = trim(slasher(inputdir))//trim(p_surf_file) if (.not.just_read) call log_param(PF, mdl, "!INPUTDIR/SURFACE_HEIGHT_IC_FILE", filename) @@ -2177,14 +2177,13 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param ! to the North/South Pole past the limits of the input data, they are extrapolated using the average ! value at the northernmost/southernmost latitude. - call horiz_interp_and_extrap_tracer(tfilename, potemp_var, 1.0, 1, & G, temp_z, mask_z, z_in, z_edges_in, missing_value_temp, reentrant_x, & - tripolar_n, homogenize, m_to_Z=US%m_to_Z, answers_2018=hor_regrid_answers_2018,ongrid=pre_gridded) + tripolar_n, homogenize, m_to_Z=US%m_to_Z, answers_2018=hor_regrid_answers_2018, ongrid=pre_gridded) call horiz_interp_and_extrap_tracer(sfilename, salin_var, 1.0, 1, & G, salt_z, mask_z, z_in, z_edges_in, missing_value_salt, reentrant_x, & - tripolar_n, homogenize, m_to_Z=US%m_to_Z, answers_2018=hor_regrid_answers_2018,ongrid=pre_gridded) + tripolar_n, homogenize, m_to_Z=US%m_to_Z, answers_2018=hor_regrid_answers_2018, ongrid=pre_gridded) kd = size(z_in,1) diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index dda892dc3e..a0f1631d6d 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -558,7 +558,7 @@ subroutine propagate_int_tide(h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, & ! Output 2-D energy loss (summed over angles) for each freq and mode do m=1,CS%NMode ; do fr=1,CS%Nfreq if (CS%id_itidal_loss_mode(fr,m) > 0 .or. CS%id_allprocesses_loss_mode(fr,m) > 0) then - itidal_loss_mode(:,:) = 0.0 ! wave-drag processes (could do others as well) + itidal_loss_mode(:,:) = 0.0 ! wave-drag processes (could do others as well) allprocesses_loss_mode(:,:) = 0.0 ! all processes summed together do a=1,CS%nAngle ; do j=js,je ; do i=is,ie itidal_loss_mode(i,j) = itidal_loss_mode(i,j) + CS%TKE_itidal_loss(i,j,a,fr,m) @@ -886,12 +886,17 @@ subroutine PPM_angular_advect(En2d, CFL_ang, Flux_En, NAngle, dt, halo_ang) !! across angles [R Z3 T-2 ~> J m-2]. ! Local variables real :: flux - real :: u_ang - real :: Angle_size - real :: I_Angle_size - real :: I_dt + real :: u_ang ! Angular propagation speed [Rad T-1 ~> Rad s-1] + real :: Angle_size ! The size of each orientation wedge in radians [Rad] + real :: I_Angle_size ! The inverse of the the orientation wedges [Rad-1] + real :: I_dt ! The inverse of the timestep [T-1 ~> s-1] + real :: aR, aL ! Left and right edge estimates of energy density [R Z3 T-2 rad-1 ~> J m-2 rad-1] + real :: dMx, dMn + real :: Ep, Ec, Em ! Mean angular energy density for three successive wedges in angular + ! orientation [R Z3 T-2 rad-1 ~> J m-2 rad-1] + real :: dA, curv_3 ! Difference and curvature of energy density [R Z3 T-2 rad-1 ~> J m-2 rad-1] + real, parameter :: oneSixth = 1.0/6.0 ! One sixth [nondim] integer :: a - real :: aR, aL, dMx, dMn, Ep, Ec, Em, dA, mA, a6 I_dt = 1 / dt Angle_size = (8.0*atan(1.0)) / (real(NAngle)) @@ -902,50 +907,55 @@ subroutine PPM_angular_advect(En2d, CFL_ang, Flux_En, NAngle, dt, halo_ang) u_ang = CFL_ang(A)*Angle_size*I_dt if (u_ang >= 0.0) then ! Implementation of PPM-H3 - Ep = En2d(a+1)*I_Angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad) - Ec = En2d(a) *I_Angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad) - Em = En2d(a-1)*I_Angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad) - aL = ( 5.*Ec + ( 2.*Em - Ep ) )/6. ! H3 estimate + ! Convert wedge-integrated energy density into angular energy densities for three successive + ! wedges around the source wedge for this flux [R Z3 T-2 rad-1 ~> J m-2 rad-1]. + Ep = En2d(a+1)*I_Angle_size + Ec = En2d(a) *I_Angle_size + Em = En2d(a-1)*I_Angle_size + ! Calculate and bound edge values of energy density. + aL = ( 5.*Ec + ( 2.*Em - Ep ) ) * oneSixth ! H3 estimate aL = max( min(Ec,Em), aL) ; aL = min( max(Ec,Em), aL) ! Bound - aR = ( 5.*Ec + ( 2.*Ep - Em ) )/6. ! H3 estimate + aR = ( 5.*Ec + ( 2.*Ep - Em ) ) * oneSixth ! H3 estimate aR = max( min(Ec,Ep), aR) ; aR = min( max(Ec,Ep), aR) ! Bound - dA = aR - aL ; mA = 0.5*( aR + aL ) + dA = aR - aL if ((Ep-Ec)*(Ec-Em) <= 0.) then - aL = Ec ; aR = Ec ! PCM for local extremum - elseif ( dA*(Ec-mA) > (dA*dA)/6. ) then - aL = 3.*Ec - 2.*aR !? - elseif ( dA*(Ec-mA) < - (dA*dA)/6. ) then - aR = 3.*Ec - 2.*aL !? + aL = Ec ; aR = Ec ! use PCM for local extremum + elseif ( 3.0*dA*(2.*Ec - (aR + aL)) > (dA*dA) ) then + aL = 3.*Ec - 2.*aR ! Flatten the profile to move the extremum to the left edge + elseif ( 3.0*dA*(2.*Ec - (aR + aL)) < - (dA*dA) ) then + aR = 3.*Ec - 2.*aL ! Flatten the profile to move the extremum to the right edge endif - a6 = 6.*Ec - 3. * (aR + aL) ! Curvature - ! CALCULATE FLUX RATE (Jm-2/s) - flux = u_ang*( aR + 0.5 * CFL_ang(A) * ( ( aL - aR ) + a6 * ( 1. - 2./3. * CFL_ang(A) ) ) ) - !flux = u_ang*( aR - 0.5 * CFL_ang(A) * ( ( aR - aL ) - a6 * ( 1. - 2./3. * CFL_ang(A) ) ) ) - ! CALCULATE AMOUNT FLUXED (Jm-2) + curv_3 = (aR + aL) - 2.0*Ec ! Curvature + ! Calculate angular flux rate [R Z3 T-3 ~> W m-2] + flux = u_ang*( aR + CFL_ang(A) * ( 0.5*(aL - aR) + curv_3 * (CFL_ang(A) - 1.5) ) ) + ! Calculate amount of energy fluxed between wedges [R Z3 T-2 ~> J m-2] Flux_En(A) = dt * flux !Flux_En(A) = (dt * I_Angle_size) * flux else ! Implementation of PPM-H3 - Ep = En2d(a+2)*I_Angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad) - Ec = En2d(a+1)*I_Angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad) - Em = En2d(a) *I_Angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad) - aL = ( 5.*Ec + ( 2.*Em - Ep ) )/6. ! H3 estimate + ! Convert wedge-integrated energy density into angular energy densities for three successive + ! wedges around the source wedge for this flux [R Z3 T-2 rad-1 ~> J m-2 rad-1]. + Ep = En2d(a+2)*I_Angle_size + Ec = En2d(a+1)*I_Angle_size + Em = En2d(a) *I_Angle_size + ! Calculate and bound edge values of energy density. + aL = ( 5.*Ec + ( 2.*Em - Ep ) ) * oneSixth ! H3 estimate aL = max( min(Ec,Em), aL) ; aL = min( max(Ec,Em), aL) ! Bound - aR = ( 5.*Ec + ( 2.*Ep - Em ) )/6. ! H3 estimate + aR = ( 5.*Ec + ( 2.*Ep - Em ) ) * oneSixth ! H3 estimate aR = max( min(Ec,Ep), aR) ; aR = min( max(Ec,Ep), aR) ! Bound - dA = aR - aL ; mA = 0.5*( aR + aL ) + dA = aR - aL if ((Ep-Ec)*(Ec-Em) <= 0.) then - aL = Ec ; aR = Ec ! PCM for local extremum - elseif ( dA*(Ec-mA) > (dA*dA)/6. ) then - aL = 3.*Ec - 2.*aR - elseif ( dA*(Ec-mA) < - (dA*dA)/6. ) then - aR = 3.*Ec - 2.*aL + aL = Ec ; aR = Ec ! use PCM for local extremum + elseif ( 3.0*dA*(2.*Ec - (aR + aL)) > (dA*dA) ) then + aL = 3.*Ec - 2.*aR ! Flatten the profile to move the extremum to the left edge + elseif ( 3.0*dA*(2.*Ec - (aR + aL)) < - (dA*dA) ) then + aR = 3.*Ec - 2.*aL ! Flatten the profile to move the extremum to the right edge endif - a6 = 6.*Ec - 3. * (aR + aL) ! Curvature - ! CALCULATE FLUX RATE (Jm-2/s) - flux = u_ang*( aR + 0.5 * CFL_ang(A) * ( ( aL - aR ) + a6 * ( 1. - 2./3. * CFL_ang(A) ) ) ) - !flux = u_ang*( aL + 0.5 * CFL_ang(A) * ( ( aR - aL ) + a6 * ( 1. - 2./3. * CFL_ang(A) ) ) ) - ! CALCULATE AMOUNT FLUXED (Jm-2) + curv_3 = (aR + aL) - 2.0*Ec ! Curvature + ! Calculate angular flux rate [R Z3 T-3 ~> W m-2] + ! Note that CFL_ang is negative here, so it looks odd compared with equivalent expressions. + flux = u_ang*( aL - CFL_ang(A) * ( 0.5*(aR - aL) + curv_3 * (-CFL_ang(A) - 1.5) ) ) + ! Calculate amount of energy fluxed between wedges [R Z3 T-2 ~> J m-2] Flux_En(A) = dt * flux !Flux_En(A) = (dt * I_Angle_size) * flux endif @@ -1014,7 +1024,7 @@ subroutine propagate(En, cn, freq, dt, G, US, CS, NAngle) ! FIND AVERAGE GROUP VELOCITY (SPEED) AT CELL CORNERS ! NOTE: THIS HAS NOT BE ADAPTED FOR REFLECTION YET (BDM)!! ! Fix indexing here later - speed(:,:) = 0 + speed(:,:) = 0.0 do J=jsh-1,jeh ; do I=ish-1,ieh f2 = G%CoriolisBu(I,J)**2 speed(I,J) = 0.25*(cn(i,j) + cn(i+1,j) + cn(i+1,j+1) + cn(i,j+1)) * & @@ -1058,21 +1068,21 @@ subroutine propagate(En, cn, freq, dt, G, US, CS, NAngle) ! Apply propagation in x-direction (reflection included) LB%jsh = jsh ; LB%jeh = jeh ; LB%ish = ish ; LB%ieh = ieh - call propagate_x(En(:,:,:), speed_x, Cgx_av(:), dCgx(:), dt, G, US, CS%nAngle, CS, LB) + call propagate_x(En, speed_x, Cgx_av, dCgx, dt, G, US, CS%nAngle, CS, LB) ! Check for energy conservation on computational domain (for debugging) - !call sum_En(G,CS,En(:,:,:),'post-propagate_x') + !call sum_En(G, CS, En, 'post-propagate_x') ! Update halos - call pass_var(En(:,:,:),G%domain) + call pass_var(En, G%domain) ! Apply propagation in y-direction (reflection included) ! LB%jsh = js ; LB%jeh = je ; LB%ish = is ; LB%ieh = ie ! Use if no teleport LB%jsh = jsh ; LB%jeh = jeh ; LB%ish = ish ; LB%ieh = ieh - call propagate_y(En(:,:,:), speed_y, Cgy_av(:), dCgy(:), dt, G, US, CS%nAngle, CS, LB) + call propagate_y(En, speed_y, Cgy_av, dCgy, dt, G, US, CS%nAngle, CS, LB) ! Check for energy conservation on computational domain (for debugging) - !call sum_En(G,CS,En(:,:,:),'post-propagate_y') + !call sum_En(G, CS, En, 'post-propagate_y') endif end subroutine propagate @@ -1084,7 +1094,7 @@ subroutine propagate_corner_spread(En, energized_wedge, NAngle, speed, dt, G, CS type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. real, dimension(G%isd:G%ied,G%jsd:G%jed), & intent(inout) :: En !< The energy density integrated over an angular - !! band [R Z3 T-2 ~> J m-2], intent in/out. + !! band [R Z3 T-2 ~> J m-2]. real, dimension(G%IsdB:G%IedB,G%Jsd:G%Jed), & intent(in) :: speed !< The magnitude of the group velocity at the cell !! corner points [L T-1 ~> m s-1]. @@ -1351,7 +1361,7 @@ subroutine propagate_x(En, speed_x, Cgx_av, dCgx, dt, G, US, Nangle, CS, LB) !! discretized wave energy spectrum. real, dimension(G%isd:G%ied,G%jsd:G%jed,Nangle), & intent(inout) :: En !< The energy density integrated over an angular - !! band [R Z3 T-2 ~> J m-2], intent in/out. + !! band [R Z3 T-2 ~> J m-2]. real, dimension(G%IsdB:G%IedB,G%jsd:G%jed), & intent(in) :: speed_x !< The magnitude of the group velocity at the !! Cu points [L T-1 ~> m s-1]. @@ -1404,18 +1414,18 @@ subroutine propagate_x(En, speed_x, Cgx_av, dCgx, dt, G, US, Nangle, CS, LB) enddo ! a-loop ! Only reflect newly arrived energy; existing energy in incident wedge is not reflected - ! and will eventually propagate out of cell. (Thid code only reflects if En > 0) - call reflect(Fdt_m(:,:,:), Nangle, CS, G, LB) - call teleport(Fdt_m(:,:,:), Nangle, CS, G, LB) - call reflect(Fdt_p(:,:,:), Nangle, CS, G, LB) - call teleport(Fdt_p(:,:,:), Nangle, CS, G, LB) + ! and will eventually propagate out of cell. (This code only reflects if En > 0.) + call reflect(Fdt_m, Nangle, CS, G, LB) + call teleport(Fdt_m, Nangle, CS, G, LB) + call reflect(Fdt_p, Nangle, CS, G, LB) + call teleport(Fdt_p, Nangle, CS, G, LB) - ! Update reflected energy (Jm-2) - do j=jsh,jeh ; do i=ish,ieh + ! Update reflected energy [R Z3 T-2 ~> J m-2] + do a=1,Nangle ; do j=jsh,jeh ; do i=ish,ieh ! if ((En(i,j,a) + G%IareaT(i,j)*(Fdt_m(i,j,a) + Fdt_p(i,j,a))) < 0.0) & ! for debugging ! call MOM_error(FATAL, "propagate_x: OutFlux>Available") - En(i,j,:) = En(i,j,:) + G%IareaT(i,j)*(Fdt_m(i,j,:) + Fdt_p(i,j,:)) - enddo ; enddo + En(i,j,a) = En(i,j,a) + G%IareaT(i,j)*(Fdt_m(i,j,a) + Fdt_p(i,j,a)) + enddo ; enddo ; enddo end subroutine propagate_x @@ -1426,7 +1436,7 @@ subroutine propagate_y(En, speed_y, Cgy_av, dCgy, dt, G, US, Nangle, CS, LB) !! discretized wave energy spectrum. real, dimension(G%isd:G%ied,G%jsd:G%jed,Nangle), & intent(inout) :: En !< The energy density integrated over an angular - !! band [R Z3 T-2 ~> J m-2], intent in/out. + !! band [R Z3 T-2 ~> J m-2]. real, dimension(G%isd:G%ied,G%JsdB:G%JedB), & intent(in) :: speed_y !< The magnitude of the group velocity at the !! Cv points [L T-1 ~> m s-1]. @@ -1486,13 +1496,13 @@ subroutine propagate_y(En, speed_y, Cgy_av, dCgy, dt, G, US, Nangle, CS, LB) enddo ! a-loop ! Only reflect newly arrived energy; existing energy in incident wedge is not reflected - ! and will eventually propagate out of cell. (Thid code only reflects if En > 0) - call reflect(Fdt_m(:,:,:), Nangle, CS, G, LB) - call teleport(Fdt_m(:,:,:), Nangle, CS, G, LB) - call reflect(Fdt_p(:,:,:), Nangle, CS, G, LB) - call teleport(Fdt_p(:,:,:), Nangle, CS, G, LB) + ! and will eventually propagate out of cell. (This code only reflects if En > 0.) + call reflect(Fdt_m, Nangle, CS, G, LB) + call teleport(Fdt_m, Nangle, CS, G, LB) + call reflect(Fdt_p, Nangle, CS, G, LB) + call teleport(Fdt_p, Nangle, CS, G, LB) - ! Update reflected energy (Jm-2) + ! Update reflected energy [R Z3 T-2 ~> J m-2] do a=1,Nangle ; do j=jsh,jeh ; do i=ish,ieh ! if ((En(i,j,a) + G%IareaT(i,j)*(Fdt_m(i,j,a) + Fdt_p(i,j,a))) < 0.0) & ! for debugging ! call MOM_error(FATAL, "propagate_y: OutFlux>Available", .true.) @@ -1521,8 +1531,7 @@ subroutine zonal_flux_En(u, h, hL, hR, uh, dt, G, US, j, ish, ieh, vol_CFL) !! the cell areas when estimating the CFL number. ! Local variables real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim]. - real :: curv_3 ! A measure of the thickness curvature over a grid length, - ! with the same units as h_in. + real :: curv_3 ! A measure of the energy density curvature over a grid length [R Z3 T-2 ~> J m-2] integer :: i do I=ish-1,ieh @@ -1566,8 +1575,7 @@ subroutine merid_flux_En(v, h, hL, hR, vh, dt, G, US, J, ish, ieh, vol_CFL) !! the CFL number. ! Local variables real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim]. - real :: curv_3 ! A measure of the thickness curvature over a grid length, - ! with the same units as h_in. + real :: curv_3 ! A measure of the energy density curvature over a grid length [R Z3 T-2 ~> J m-2] integer :: i do i=ish,ieh @@ -1603,18 +1611,18 @@ subroutine reflect(En, NAngle, CS, G, LB) type(loop_bounds_type), intent(in) :: LB !< A structure with the active energy loop bounds. ! Local variables real, dimension(G%isd:G%ied,G%jsd:G%jed) :: angle_c - ! angle of boudary wrt equator + ! angle of boundary wrt equator [rad] real, dimension(G%isd:G%ied,G%jsd:G%jed) :: part_refl ! fraction of wave energy reflected - ! values should collocate with angle_c + ! values should collocate with angle_c [nondim] logical, dimension(G%isd:G%ied,G%jsd:G%jed) :: ridge ! tags of cells with double reflection - real :: TwoPi ! 2*pi - real :: Angle_size ! size of beam wedge (rad) - real :: angle_wall ! angle of coast/ridge/shelf wrt equator - real, dimension(1:NAngle) :: angle_i ! angle of incident ray wrt equator - real :: angle_r ! angle of reflected ray wrt equator + real :: TwoPi ! 2*pi = 6.2831853... [nondim] + real :: Angle_size ! size of beam wedge [rad] + real :: angle_wall ! angle of coast/ridge/shelf wrt equator [rad] + real, dimension(1:NAngle) :: angle_i ! angle of incident ray wrt equator [rad] + real :: angle_r ! angle of reflected ray wrt equator [rad] real, dimension(1:Nangle) :: En_reflected integer :: i, j, a, a_r, na !integer :: isd, ied, jsd, jed ! start and end local indices on data domain @@ -1623,7 +1631,6 @@ subroutine reflect(En, NAngle, CS, G, LB) ! (values exclude halos) integer :: ish, ieh, jsh, jeh ! start and end local indices on data domain ! leaving out outdated halo points (march in) - integer :: id_g, jd_g ! global (decomp-invar) indices !isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec @@ -1643,59 +1650,54 @@ subroutine reflect(En, NAngle, CS, G, LB) ridge = CS%refl_dbl En_reflected(:) = 0.0 - !do j=jsc-1,jec+1 - do j=jsh,jeh - !do i=isc-1,iec+1 - do i=ish,ieh - ! jd_g = j + G%jdg_offset ; id_g = i + G%idg_offset - ! redistribute energy in angular space if ray will hit boundary - ! i.e., if energy is in a reflecting cell - if (angle_c(i,j) /= CS%nullangle) then - do a=1,NAngle - if (En(i,j,a) > 0.0) then - ! if ray is incident, keep specified boundary angle - if (sin(angle_i(a) - angle_c(i,j)) >= 0.0) then - angle_wall = angle_c(i,j) - ! if ray is not incident but in ridge cell, use complementary angle - elseif (ridge(i,j)) then - angle_wall = angle_c(i,j) + 0.5*TwoPi - if (angle_wall > TwoPi) then - angle_wall = angle_wall - TwoPi*floor(abs(angle_wall)/TwoPi) - elseif (angle_wall < 0.0) then - angle_wall = angle_wall + TwoPi*ceiling(abs(angle_wall)/TwoPi) - endif - ! if ray is not incident and not in a ridge cell, keep specified angle - else - angle_wall = angle_c(i,j) - endif - ! do reflection - if (sin(angle_i(a) - angle_wall) >= 0.0) then - angle_r = 2.0*angle_wall - angle_i(a) - if (angle_r > TwoPi) then - angle_r = angle_r - TwoPi*floor(abs(angle_r)/TwoPi) - elseif (angle_r < 0.0) then - angle_r = angle_r + TwoPi*ceiling(abs(angle_r)/TwoPi) - endif - a_r = nint(angle_r/Angle_size) + 1 - do while (a_r > Nangle) ; a_r = a_r - Nangle ; enddo - if (a /= a_r) then - En_reflected(a_r) = part_refl(i,j)*En(i,j,a) - En(i,j,a) = (1.0-part_refl(i,j))*En(i,j,a) - endif - endif + do j=jsh,jeh ; do i=ish,ieh + ! redistribute energy in angular space if ray will hit boundary + ! i.e., if energy is in a reflecting cell + if (angle_c(i,j) /= CS%nullangle) then + do a=1,NAngle ; if (En(i,j,a) > 0.0) then + if (sin(angle_i(a) - angle_c(i,j)) >= 0.0) then + ! if ray is incident, keep specified boundary angle + angle_wall = angle_c(i,j) + elseif (ridge(i,j)) then + ! if ray is not incident but in ridge cell, use complementary angle + angle_wall = angle_c(i,j) + 0.5*TwoPi + if (angle_wall > TwoPi) then + angle_wall = angle_wall - TwoPi*floor(abs(angle_wall)/TwoPi) + elseif (angle_wall < 0.0) then + angle_wall = angle_wall + TwoPi*ceiling(abs(angle_wall)/TwoPi) endif - enddo ! a-loop - En(i,j,:) = En(i,j,:) + En_reflected(:) - En_reflected(:) = 0.0 - endif - enddo ! i-loop - enddo ! j-loop + else + ! if ray is not incident and not in a ridge cell, keep specified angle + angle_wall = angle_c(i,j) + endif + + ! do reflection + if (sin(angle_i(a) - angle_wall) >= 0.0) then + angle_r = 2.0*angle_wall - angle_i(a) + if (angle_r > TwoPi) then + angle_r = angle_r - TwoPi*floor(abs(angle_r)/TwoPi) + elseif (angle_r < 0.0) then + angle_r = angle_r + TwoPi*ceiling(abs(angle_r)/TwoPi) + endif + a_r = nint(angle_r/Angle_size) + 1 + do while (a_r > Nangle) ; a_r = a_r - Nangle ; enddo + if (a /= a_r) then + En_reflected(a_r) = part_refl(i,j)*En(i,j,a) + En(i,j,a) = (1.0-part_refl(i,j))*En(i,j,a) + endif + endif + endif ; enddo ! a-loop + do a=1,NAngle + En(i,j,a) = En(i,j,a) + En_reflected(a) + En_reflected(a) = 0.0 + enddo ! a-loop + endif + enddo ; enddo ! i- and j-loops ! Check to make sure no energy gets onto land (only run for debugging) ! do a=1,NAngle ; do j=jsc,jec ; do i=isc,iec ! if (En(i,j,a) > 0.001 .and. G%mask2dT(i,j) == 0) then - ! jd_g = j + G%jdg_offset ; id_g = i + G%idg_offset - ! write (mesg,*) 'En=', En(i,j,a), 'a=', a, 'ig_g=',id_g, 'jg_g=',jd_g + ! write (mesg,*) 'En=', En(i,j,a), 'a=', a, 'ig_g=',i+G%idg_offset, 'jg_g=',j+G%jdg_offset ! call MOM_error(FATAL, "reflect: Energy detected out of bounds: "//trim(mesg), .true.) ! endif ! enddo ; enddo ; enddo @@ -1717,17 +1719,17 @@ subroutine teleport(En, NAngle, CS, G, LB) type(loop_bounds_type), intent(in) :: LB !< A structure with the active energy loop bounds. ! Local variables real, dimension(G%isd:G%ied,G%jsd:G%jed) :: angle_c - ! angle of boudary wrt equator + ! angle of boundary wrt equator [rad] real, dimension(G%isd:G%ied,G%jsd:G%jed) :: part_refl ! fraction of wave energy reflected - ! values should collocate with angle_c + ! values should collocate with angle_c [nondim] logical, dimension(G%isd:G%ied,G%jsd:G%jed) :: pref_cell ! flag for partial reflection logical, dimension(G%isd:G%ied,G%jsd:G%jed) :: ridge - ! tags of cells with double reflection - real :: TwoPi ! size of beam wedge (rad) - real :: Angle_size ! size of beam wedge (rad) - real, dimension(1:NAngle) :: angle_i ! angle of incident ray wrt equator + ! tags of cells with double reflection + real :: TwoPi ! 2*pi = 6.2831853... [nondim] + real :: Angle_size ! size of beam wedge [rad] + real, dimension(1:NAngle) :: angle_i ! angle of incident ray wrt equator [rad] real, dimension(1:NAngle) :: cos_angle, sin_angle real :: En_tele ! energy to be "teleported" [R Z3 T-2 ~> J m-2] character(len=160) :: mesg ! The text of an error message @@ -2295,8 +2297,8 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) CS%TKE_itidal_loss(:,:,:,:,:) = 0.0 allocate(CS%TKE_Froude_loss(isd:ied,jsd:jed,num_angle,num_freq,num_mode)) CS%TKE_Froude_loss(:,:,:,:,:) = 0.0 - allocate(CS%tot_leak_loss(isd:ied,jsd:jed)) ; CS%tot_leak_loss(:,:) = 0.0 - allocate(CS%tot_quad_loss(isd:ied,jsd:jed) ) ; CS%tot_quad_loss(:,:) = 0.0 + allocate(CS%tot_leak_loss(isd:ied,jsd:jed)) ; CS%tot_leak_loss(:,:) = 0.0 + allocate(CS%tot_quad_loss(isd:ied,jsd:jed) ) ; CS%tot_quad_loss(:,:) = 0.0 allocate(CS%tot_itidal_loss(isd:ied,jsd:jed)) ; CS%tot_itidal_loss(:,:) = 0.0 allocate(CS%tot_Froude_loss(isd:ied,jsd:jed)) ; CS%tot_Froude_loss(:,:) = 0.0 diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index b791535ed1..438236fc2a 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -171,7 +171,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ integer :: i, j, k, col, total_sponge_cols, total_sponge_cols_u, total_sponge_cols_v character(len=10) :: remapScheme if (associated(CS)) then - call MOM_error(WARNING, "initialize_sponge called with an associated "// & + call MOM_error(WARNING, "initialize_ALE_sponge_fixed called with an associated "// & "control structure.") return endif @@ -260,14 +260,14 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ if (CS%sponge_uv) then - allocate(data_hu(G%isdB:G%iedB,G%jsd:G%jed,nz_data)); data_hu(:,:,:)=0.0 - allocate(data_hv(G%isd:G%ied,G%jsdB:G%jedB,nz_data)); data_hv(:,:,:)=0.0 - allocate(Iresttime_u(G%isdB:G%iedB,G%jsd:G%jed)); Iresttime_u(:,:)=0.0 - allocate(Iresttime_v(G%isd:G%ied,G%jsdB:G%jedB)); Iresttime_v(:,:)=0.0 + allocate(data_hu(G%isdB:G%iedB,G%jsd:G%jed,nz_data)) ; data_hu(:,:,:) = 0.0 + allocate(data_hv(G%isd:G%ied,G%jsdB:G%jedB,nz_data)) ; data_hv(:,:,:) = 0.0 + allocate(Iresttime_u(G%isdB:G%iedB,G%jsd:G%jed)) ; Iresttime_u(:,:) = 0.0 + allocate(Iresttime_v(G%isd:G%ied,G%jsdB:G%jedB)) ; Iresttime_v(:,:) = 0.0 ! u points CS%num_col_u = 0 ; !CS%fldno_u = 0 - do j=CS%jsc,CS%jec; do I=CS%iscB,CS%iecB + do j=CS%jsc,CS%jec ; do I=CS%iscB,CS%iecB data_hu(I,j,:) = 0.5 * (data_h(i,j,:) + data_h(i+1,j,:)) Iresttime_u(I,j) = 0.5 * (Iresttime(i,j) + Iresttime(i+1,j)) if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) & @@ -276,9 +276,9 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ if (CS%num_col_u > 0) then - allocate(CS%Iresttime_col_u(CS%num_col_u)) ; CS%Iresttime_col_u = 0.0 - allocate(CS%col_i_u(CS%num_col_u)) ; CS%col_i_u = 0 - allocate(CS%col_j_u(CS%num_col_u)) ; CS%col_j_u = 0 + allocate(CS%Iresttime_col_u(CS%num_col_u)) ; CS%Iresttime_col_u(:) = 0.0 + allocate(CS%col_i_u(CS%num_col_u)) ; CS%col_i_u(:) = 0 + allocate(CS%col_j_u(CS%num_col_u)) ; CS%col_j_u(:) = 0 ! pass indices, restoring time to the CS structure col = 1 @@ -286,7 +286,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) then CS%col_i_u(col) = i ; CS%col_j_u(col) = j CS%Iresttime_col_u(col) = Iresttime_u(i,j) - col = col +1 + col = col + 1 endif enddo ; enddo @@ -323,7 +323,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ if ((Iresttime_v(i,J)>0.0) .and. (G%mask2dCv(i,J)>0)) then CS%col_i_v(col) = i ; CS%col_j_v(col) = j CS%Iresttime_col_v(col) = Iresttime_v(i,j) - col = col +1 + col = col + 1 endif enddo ; enddo @@ -415,7 +415,7 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, param_file, CS) character(len=10) :: remapScheme if (associated(CS)) then - call MOM_error(WARNING, "initialize_sponge called with an associated "// & + call MOM_error(WARNING, "initialize_ALE_sponge_varying called with an associated "// & "control structure.") return endif @@ -486,8 +486,8 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, param_file, CS) call log_param(param_file, mdl, "!Total sponge columns at h points", total_sponge_cols, & "The total number of columns where sponges are applied at h points.") if (CS%sponge_uv) then - allocate(Iresttime_u(G%isdB:G%iedB,G%jsd:G%jed)); Iresttime_u(:,:)=0.0 - allocate(Iresttime_v(G%isd:G%ied,G%jsdB:G%jedB)); Iresttime_v(:,:)=0.0 + allocate(Iresttime_u(G%isdB:G%iedB,G%jsd:G%jed)) ; Iresttime_u(:,:) = 0.0 + allocate(Iresttime_v(G%isd:G%ied,G%jsdB:G%jedB)) ; Iresttime_v(:,:) = 0.0 ! u points CS%num_col_u = 0 ; !CS%fldno_u = 0 do j=CS%jsc,CS%jec; do I=CS%iscB,CS%iecB @@ -578,7 +578,7 @@ subroutine set_up_ALE_sponge_field_fixed(sp_val, G, f_ptr, CS) if (CS%fldno > MAX_FIELDS_) then write(mesg,'("Increase MAX_FIELDS_ to at least ",I3," in MOM_memory.h or decrease & &the number of fields to be damped in the call to & - &initialize_sponge." )') CS%fldno + &initialize_ALE_sponge." )') CS%fldno call MOM_error(FATAL,"set_up_ALE_sponge_field: "//mesg) endif @@ -605,8 +605,8 @@ subroutine set_up_ALE_sponge_field_varying(filename, fieldname, Time, G, GV, US, type(time_type), intent(in) :: Time !< The current model time type(ocean_grid_type), intent(in) :: G !< Grid structure (in). type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & target, intent(in) :: f_ptr !< Pointer to the field to be damped (in). type(ALE_sponge_CS), pointer :: CS !< Sponge control structure (in/out). @@ -634,7 +634,7 @@ subroutine set_up_ALE_sponge_field_varying(filename, fieldname, Time, G, GV, US, if (CS%fldno > MAX_FIELDS_) then write(mesg,'("Increase MAX_FIELDS_ to at least ",I3," in MOM_memory.h or decrease & &the number of fields to be damped in the call to & - &initialize_sponge." )') CS%fldno + &initialize_ALE_sponge." )') CS%fldno call MOM_error(FATAL,"set_up_ALE_sponge_field: "//mesg) endif ! get a unique time interp id for this field. If sponge data is ongrid, then setup @@ -788,11 +788,11 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure (in). type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(inout) :: h !< Layer thickness [H ~> m or kg m-2] (in) real, intent(in) :: dt !< The amount of time covered by this call [T ~> s]. type(ALE_sponge_CS), pointer :: CS !< A pointer to the control structure for this module - !! that is set by a previous call to initialize_sponge (in). + !! that is set by a previous call to initialize_ALE_sponge (in). type(time_type), optional, intent(in) :: Time !< The current model date real :: damp ! The timestep times the local damping coefficient [nondim]. @@ -833,8 +833,8 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) nz_data = CS%Ref_val(m)%nz_data allocate(sp_val(G%isd:G%ied,G%jsd:G%jed,1:nz_data)) allocate(mask_z(G%isd:G%ied,G%jsd:G%jed,1:nz_data)) - sp_val(:,:,:)=0.0 - mask_z(:,:,:)=0.0 + sp_val(:,:,:) = 0.0 + mask_z(:,:,:) = 0.0 call horiz_interp_and_extrap_tracer(CS%Ref_val(m)%id, Time, 1.0, G, sp_val, mask_z, z_in, & z_edges_in, missing_value, .true., .false., .false., & spongeOnGrid=CS%SpongeDataOngrid, m_to_Z=US%m_to_Z, & @@ -1003,17 +1003,20 @@ end subroutine apply_ALE_sponge !> Rotate the ALE sponge fields from the input to the model index map. subroutine rotate_ALE_sponge(sponge_in, G_in, sponge, G, turns, param_file) - type(ALE_sponge_CS), intent(in) :: sponge_in - type(ocean_grid_type), intent(in) :: G_in - type(ALE_sponge_CS), pointer :: sponge - type(ocean_grid_type), intent(in) :: G - integer, intent(in) :: turns - type(param_file_type), intent(in) :: param_file + type(ALE_sponge_CS), intent(in) :: sponge_in !< The control structure for this module with the + !! original grid rotation + type(ocean_grid_type), intent(in) :: G_in !< The ocean's grid structure with the original rotation. + type(ALE_sponge_CS), pointer :: sponge !< A pointer to the control that will be set up with + !! the new grid rotation + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure with the new rotation. + integer, intent(in) :: turns !< The number of 90-degree turns between grids + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file + !! to parse for model parameter values. ! First part: Index construction ! 1. Reconstruct Iresttime(:,:) from sponge_in ! 2. rotate Iresttime(:,:) - ! 3. Call initialize_sponge using new grid and rotated Iresttime(:,:) + ! 3. Call initialize_ALE_sponge using new grid and rotated Iresttime(:,:) ! All the index adjustment should follow from the Iresttime rotation real, dimension(:,:), allocatable :: Iresttime_in, Iresttime @@ -1040,7 +1043,7 @@ subroutine rotate_ALE_sponge(sponge_in, G_in, sponge, G, turns, param_file) endif ! Re-populate the 2D Iresttime and data_h arrays on the original grid - do c = 1, sponge_in%num_col + do c=1,sponge_in%num_col c_i = sponge_in%col_i(c) c_j = sponge_in%col_j(c) Iresttime_in(c_i, c_j) = sponge_in%Iresttime_col(c) @@ -1076,24 +1079,25 @@ subroutine rotate_ALE_sponge(sponge_in, G_in, sponge, G, turns, param_file) allocate(sp_val(G%isd:G%ied, G%jsd:G%jed, nz_data)) endif - do n = 1, sponge_in%fldno + do n=1,sponge_in%fldno ! Assume that tracers are pointers and are remapped in other functions(?) sp_ptr => sponge_in%var(n)%p - sp_val_in(:,:,:) = 0.0 - do c = 1, sponge_in%num_col - c_i = sponge_in%col_i(c) - c_j = sponge_in%col_j(c) - if (fixed_sponge) then + if (fixed_sponge) then + sp_val_in(:,:,:) = 0.0 + do c = 1, sponge_in%num_col + c_i = sponge_in%col_i(c) + c_j = sponge_in%col_j(c) do k = 1, nz_data sp_val_in(c_i, c_j, k) = sponge_in%Ref_val(n)%p(k,c) enddo - endif - enddo + enddo + + call rotate_array(sp_val_in, turns, sp_val) - call rotate_array(sp_val_in, turns, sp_val) - if (fixed_sponge) then ! NOTE: This points sp_val with the unrotated field. See note below. call set_up_ALE_sponge_field(sp_val, G, sp_ptr, sponge) + + deallocate(sp_val_in) else ! We don't want to repeat FMS init in set_up_ALE_sponge_field_varying() ! (time_interp_external_init, init_external_field, etc), so we manually @@ -1123,11 +1127,6 @@ subroutine rotate_ALE_sponge(sponge_in, G_in, sponge, G, turns, param_file) endif enddo - if (fixed_sponge) then - deallocate(sp_val_in) - deallocate(sp_val) - endif - ! TODO: var_u and var_v sponge dampling is not yet supported. if (associated(sponge_in%var_u%p) .or. associated(sponge_in%var_v%p)) & call MOM_error(FATAL, "Rotation of ALE sponge velocities is not yet " & @@ -1144,17 +1143,22 @@ end subroutine rotate_ALE_sponge ! TODO: This function solely exists to replace field pointers in the sponge ! after rotation. This function is part of a temporary solution until ! something more robust is developed. -subroutine update_ALE_sponge_field(sponge, p_old, p_new) - type(ALE_sponge_CS), pointer :: sponge - real, dimension(:,:,:), target, intent(in) :: p_old - real, dimension(:,:,:), target, intent(in) :: p_new +subroutine update_ALE_sponge_field(sponge, p_old, G, GV, p_new) + type(ALE_sponge_CS), pointer :: sponge !< A pointer to the control structure for this module + !! that is set by a previous call to initialize_ALE_sponge. + real, dimension(:,:,:), & + target, intent(in) :: p_old !< The previous array of target values + type(ocean_grid_type), intent(in) :: G !< The updated ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + target, intent(in) :: p_new !< The new array of target values integer :: n - do n = 1, sponge%fldno - if (associated(sponge%var(n)%p, p_old)) & - sponge%var(n)%p => p_new + do n=1,sponge%fldno + if (associated(sponge%var(n)%p, p_old)) sponge%var(n)%p => p_new enddo + end subroutine update_ALE_sponge_field @@ -1163,7 +1167,7 @@ end subroutine update_ALE_sponge_field !> This subroutine deallocates any memory associated with the ALE_sponge module. subroutine ALE_sponge_end(CS) type(ALE_sponge_CS), pointer :: CS !< A pointer to the control structure that is - !! set by a previous call to initialize_sponge. + !! set by a previous call to initialize_ALE_sponge. integer :: m diff --git a/src/parameterizations/vertical/MOM_CVMix_shear.F90 b/src/parameterizations/vertical/MOM_CVMix_shear.F90 index 6f1a629ab4..f099305f0c 100644 --- a/src/parameterizations/vertical/MOM_CVMix_shear.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_shear.F90 @@ -88,10 +88,11 @@ subroutine calculate_CVMix_shear(u_H, v_H, h, tv, kd, kv, G, GV, US, CS ) real, dimension(G%ke+1) :: Ri_Grad !< Gradient Richardson number [nondim] real, dimension(G%ke+1) :: Kvisc !< Vertical viscosity at interfaces [m2 s-1] real, dimension(G%ke+1) :: Kdiff !< Diapycnal diffusivity at interfaces [m2 s-1] - real, parameter :: epsln = 1.e-10 !< Threshold to identify vanished layers [m] + real :: epsln !< Threshold to identify vanished layers [H ~> m or kg m-2] ! some constants GoRho = US%L_to_Z**2 * GV%g_Earth / GV%Rho0 + epsln = 1.e-10 * GV%m_to_H do j = G%jsc, G%jec do i = G%isc, G%iec @@ -148,7 +149,6 @@ subroutine calculate_CVMix_shear(u_H, v_H, h, tv, kd, kv, G, GV, US, CS ) if (CS%smooth_ri) then ! 1) fill Ri_grad in vanished layers with adjacent value - !### For dimensional consistency, epsln needs to be epsln*GV%m_to_H. do k = 2, G%ke if (h(i,j,k) <= epsln) Ri_grad(k) = Ri_grad(k-1) enddo diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index d753afc97b..0bd3138670 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -2273,11 +2273,8 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e endif - ! This block sets ea, eb from Kd or Kd_int. - ! Otherwise, call entrainment_diffusive() which sets ea and eb - ! based on KD and target densities (ie. does remapping as well). - ! When not using ALE, calculate layer entrainments/detrainments from - ! diffusivities and differences between layer and target densities + ! Calculate layer entrainments and detrainments from diffusivities and differences between + ! layer and target densities (i.e. do remapping as well as diffusion). call cpu_clock_begin(id_clock_entrain) ! Calculate appropriately limited diapycnal mass fluxes to account ! for diapycnal diffusion and advection. Sets: ea, eb. Changes: kb @@ -3688,7 +3685,8 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di ! False. CS%use_CVMix_conv = CVMix_conv_init(Time, G, GV, US, param_file, diag, CS%CVMix_conv_csp) - call entrain_diffusive_init(Time, G, GV, US, param_file, diag, CS%entrain_diffusive_CSp) + call entrain_diffusive_init(Time, G, GV, US, param_file, diag, CS%entrain_diffusive_CSp, & + just_read_params=CS%useALEalgorithm) ! initialize the geothermal heating module if (CS%use_geothermal) & @@ -3763,12 +3761,10 @@ subroutine diabatic_driver_end(CS) call entrain_diffusive_end(CS%entrain_diffusive_CSp) call set_diffusivity_end(CS%set_diff_CSp) - if (CS%useKPP) then + if (CS%useKPP) then deallocate( CS%KPP_buoy_flux ) deallocate( CS%KPP_temp_flux ) deallocate( CS%KPP_salt_flux ) - endif - if (CS%useKPP) then deallocate( CS%KPP_NLTheat ) deallocate( CS%KPP_NLTscalar ) call KPP_end(CS%KPP_CSp) diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index a9e68736e7..25e1f80ff0 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -4,6 +4,7 @@ module MOM_energetic_PBL ! This file is part of MOM6. See LICENSE.md for the license. use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE +use MOM_coms, only : EFP_type, real_to_EFP, EFP_to_real, operator(+), assignment(=), EFP_sum_across_PEs use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_alloc use MOM_diag_mediator, only : time_type, diag_ctrl use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type @@ -33,12 +34,12 @@ module MOM_energetic_PBL type, public :: energetic_PBL_CS ; private !/ Constants - real :: VonKar = 0.41 !< The von Karman coefficient. This should be runtime, but because - !! it is runtime in KPP and set to 0.4 it might change answers. - real :: omega !< The Earth's rotation rate [T-1 ~> s-1]. - real :: omega_frac !< When setting the decay scale for turbulence, use this fraction of - !! the absolute rotation rate blended with the local value of f, as - !! sqrt((1-of)*f^2 + of*4*omega^2) [nondim]. + real :: VonKar = 0.41 !< The von Karman coefficient. This should be a runtime parameter, + !! but because it is set to 0.4 at runtime in KPP it might change answers. + real :: omega !< The Earth's rotation rate [T-1 ~> s-1]. + real :: omega_frac !< When setting the decay scale for turbulence, use this fraction of + !! the absolute rotation rate blended with the local value of f, as + !! sqrt((1-omega_frac)*f^2 + omega_frac*4*omega^2) [nondim]. !/ Convection related terms real :: nstar !< The fraction of the TKE input to the mixed layer available to drive @@ -47,9 +48,14 @@ module MOM_energetic_PBL !! TKE produced by buoyancy. !/ Mixing Length terms - logical :: Use_MLD_iteration=.false. !< False to use old ePBL method. - logical :: MLD_iteration_guess=.false. !< False to default to guessing half the - !! ocean depth for the iteration. + logical :: Use_MLD_iteration !< If true, use the proximity to the bottom of the actively turbulent + !! surface boundary layer to constrain the mixing lengths. + logical :: MLD_iteration_guess !< False to default to guessing half the + !! ocean depth for the first iteration. + logical :: MLD_bisection !< If true, use bisection with the iterative determination of the + !! self-consistent mixed layer depth. Otherwise use the false position + !! after a maximum and minimum bound have been evaluated and the + !! returned value from the previous guess or bisection before this. integer :: max_MLD_its !< The maximum number of iterations that can be used to find a !! self-consistent mixed layer depth with Use_MLD_iteration. real :: MixLenExponent !< Exponent in the mixing length shape-function. @@ -179,6 +185,8 @@ module MOM_energetic_PBL LA, & !< Langmuir number [nondim] LA_MOD !< Modified Langmuir number [nondim] + type(EFP_type), dimension(2) :: sum_its !< The total number of iterations and columns worked on + real, allocatable, dimension(:,:,:) :: & Velocity_Scale, & !< The velocity scale used in getting Kd [Z T-1 ~> m s-1] Mixing_Length !< The length scale used in getting Kd [Z ~> m] @@ -215,6 +223,8 @@ module MOM_energetic_PBL character*(20), parameter :: ADDITIVE_STRING = "ADDITIVE" !>@} +logical :: report_avg_its = .false. !< Report the average number of ePBL iterations for debugging. + !> A type for conveniently passing around ePBL diagnostics for a column. type, public :: ePBL_column_diags ; private !>@{ Local column copies of energy change diagnostics, all in [R Z3 T-3 ~> W m-2]. @@ -755,6 +765,8 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs ! manner giving a usable guess. When it does fail, it is due to convection ! within the boundary layer. Likely, a new method e.g. surface_disconnect, ! can improve this. + real :: dMLD_min ! The change in diagnosed mixed layer depth when the guess is min_MLD [Z ~> m] + real :: dMLD_max ! The change in diagnosed mixed layer depth when the guess is max_MLD [Z ~> m] logical :: FIRST_OBL ! Flag for computing "found" Mixing layer depth logical :: OBL_converged ! Flag for convergence of MLD integer :: OBL_it ! Iteration counter @@ -827,8 +839,10 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs !/The following lines are for the iteration over MLD ! max_MLD will initialized as ocean bottom depth max_MLD = 0.0 ; do k=1,nz ; max_MLD = max_MLD + h(k)*GV%H_to_Z ; enddo - !min_MLD will initialize as 0. + ! min_MLD will be initialized to 0. min_MLD = 0.0 + ! Set values of the wrong signs to indicate that these changes are not based on valid estimates + dMLD_min = -1.0*US%m_to_Z ; dMLD_max = 1.0*US%m_to_Z ! If no first guess is provided for MLD, try the middle of the water column if (MLD_guess <= min_MLD) MLD_guess = 0.5 * (min_MLD + max_MLD) @@ -1408,17 +1422,37 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs !New method uses ML_DEPTH as computed in ePBL routine MLD_found = MLD_output - if (MLD_found - CS%MLD_tol > MLD_guess) then - min_MLD = MLD_guess - elseif (abs(MLD_guess - MLD_found) < CS%MLD_tol) then + if (MLD_found - MLD_guess > CS%MLD_tol) then + min_MLD = MLD_guess ; dMLD_min = MLD_found - MLD_guess + elseif (abs(MLD_found - MLD_guess) < CS%MLD_tol) then OBL_converged = .true. ! Break convergence loop - else - max_MLD = MLD_guess ! We know this guess was too deep + else ! We know this guess was too deep + max_MLD = MLD_guess ; dMLD_max = MLD_found - MLD_guess ! < -CS%MLD_tol endif - ! For next pass, guess average of minimum and maximum values. - !### We should try using the false position method instead of simple bisection. - MLD_guess = 0.5*(min_MLD + max_MLD) + if (.not.OBL_converged) then ; if (CS%MLD_bisection) then + ! For the next pass, guess the average of the minimum and maximum values. + MLD_guess = 0.5*(min_MLD + max_MLD) + else ! Try using the false position method or the returned value instead of simple bisection. + ! Taking the occasional step with MLD_output empirically helps to converge faster. + if ((dMLD_min > 0.0) .and. (dMLD_max < 0.0) .and. (OBL_it > 2) .and. (mod(OBL_it-1,4)>0)) then + ! Both bounds have valid change estimates and are probably in the range of possible outputs. + MLD_Guess = (dMLD_min*max_MLD - dMLD_max*min_MLD) / (dMLD_min - dMLD_max) + elseif ((MLD_found > min_MLD) .and. (MLD_found < max_MLD)) then + ! The output MLD_found is an interesting guess, as it likely to bracket the true solution + ! along with the previous value of MLD_guess and to be close to the solution. + MLD_guess = MLD_found + else ! Bisect if the other guesses would be out-of-bounds. This does not happen much. + MLD_guess = 0.5*(min_MLD + max_MLD) + endif + endif ; endif + endif + if ((OBL_converged) .or. (OBL_it==CS%Max_MLD_Its)) then + if (report_avg_its) then + CS%sum_its(1) = CS%sum_its(1) + real_to_EFP(real(OBL_it)) + CS%sum_its(2) = CS%sum_its(2) + real_to_EFP(1.0) + endif + exit endif enddo ! Iteration loop for converged boundary layer thickness. if (CS%Use_LT) then @@ -2131,16 +2165,22 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) endif call get_param(param_file, mdl, "MLD_ITERATION_GUESS", CS%MLD_ITERATION_GUESS, & - "A logical that specifies whether or not to use the "//& - "previous timestep MLD as a first guess in the MLD iteration. "//& - "The default is false to facilitate reproducibility.", default=.false.) + "If true, use the previous timestep MLD as a first guess in the MLD iteration, "//& + "otherwise use half the ocean depth as the first guess of the boundary layer "//& + "depth. The default is false to facilitate reproducibility.", & + default=.false., do_not_log=.not.CS%Use_MLD_iteration) call get_param(param_file, mdl, "EPBL_MLD_TOLERANCE", CS%MLD_tol, & "The tolerance for the iteratively determined mixed "//& "layer depth. This is only used with USE_MLD_ITERATION.", & - units="meter", default=1.0, scale=US%m_to_Z) + units="meter", default=1.0, scale=US%m_to_Z, do_not_log=.not.CS%Use_MLD_iteration) + call get_param(param_file, mdl, "EPBL_MLD_BISECTION", CS%MLD_bisection, & + "If true, use bisection with the iterative determination of the self-consistent "//& + "mixed layer depth. Otherwise use the false position after a maximum and minimum "//& + "bound have been evaluated and the returned value or bisection before this.", & + default=.true., do_not_log=.not.CS%Use_MLD_iteration) !### The default should become false. call get_param(param_file, mdl, "EPBL_MLD_MAX_ITS", CS%max_MLD_its, & "The maximum number of iterations that can be used to find a self-consistent "//& - "mixed layer depth. For now, due to the use of bisection, the maximum number "//& + "mixed layer depth. If EPBL_MLD_BISECTION is true, the maximum number "//& "iteractions needed is set by Depth/2^MAX_ITS < EPBL_MLD_TOLERANCE.", & default=20, do_not_log=.not.CS%Use_MLD_iteration) if (.not.CS%Use_MLD_iteration) CS%Max_MLD_Its = 1 @@ -2339,6 +2379,10 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) "If true, temperature and salinity are used as state "//& "variables.", default=.true.) + if (report_avg_its) then + CS%sum_its(1) = real_to_EFP(0.0) ; CS%sum_its(2) = real_to_EFP(0.0) + endif + if (max(CS%id_TKE_wind, CS%id_TKE_MKE, CS%id_TKE_conv, & CS%id_TKE_mixing, CS%id_TKE_mech_decay, CS%id_TKE_forcing, & CS%id_TKE_conv_decay) > 0) then @@ -2370,6 +2414,9 @@ subroutine energetic_PBL_end(CS) type(energetic_PBL_CS), pointer :: CS !< Energetic_PBL control structure that !! will be deallocated in this subroutine. + character(len=256) :: mesg + real :: avg_its + if (.not.associated(CS)) return if (allocated(CS%ML_depth)) deallocate(CS%ML_depth) @@ -2387,6 +2434,14 @@ subroutine energetic_PBL_end(CS) if (allocated(CS%Mixing_Length)) deallocate(CS%Mixing_Length) if (allocated(CS%Velocity_Scale)) deallocate(CS%Velocity_Scale) + if (report_avg_its) then + call EFP_sum_across_PEs(CS%sum_its, 2) + + avg_its = EFP_to_real(CS%sum_its(1)) / EFP_to_real(CS%sum_its(2)) + write (mesg,*) "Average ePBL iterations = ", avg_its + call MOM_mesg(mesg) + endif + deallocate(CS) end subroutine energetic_PBL_end diff --git a/src/parameterizations/vertical/MOM_entrain_diffusive.F90 b/src/parameterizations/vertical/MOM_entrain_diffusive.F90 index 4e30756f7b..1be3421534 100644 --- a/src/parameterizations/vertical/MOM_entrain_diffusive.F90 +++ b/src/parameterizations/vertical/MOM_entrain_diffusive.F90 @@ -2081,7 +2081,7 @@ end subroutine find_maxF_kb !> This subroutine initializes the parameters and memory associated with the !! entrain_diffusive module. -subroutine entrain_diffusive_init(Time, G, GV, US, param_file, diag, CS) +subroutine entrain_diffusive_init(Time, G, GV, US, param_file, diag, CS, just_read_params) type(time_type), intent(in) :: Time !< The current model time. type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. @@ -2092,18 +2092,15 @@ subroutine entrain_diffusive_init(Time, G, GV, US, param_file, diag, CS) !! output. type(entrain_diffusive_CS), pointer :: CS !< A pointer that is set to point to the control !! structure. -! for this module -! Arguments: Time - The current model time. -! (in) G - The ocean's grid structure. -! (in) GV - The ocean's vertical grid structure. -! (in) param_file - A structure indicating the open file to parse for -! model parameter values. -! (in) diag - A structure that is used to regulate diagnostic output. -! (in/out) CS - A pointer that is set to point to the control structure -! for this module + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters logging them or registering + !! any diagnostics + + ! Local variables real :: decay_length, dt, Kd -! This include declares and sets the variable "version". -#include "version_variable.h" + logical :: just_read ! If true, just read parameters but do nothing else. + ! This include declares and sets the variable "version". +# include "version_variable.h" character(len=40) :: mdl = "MOM_entrain_diffusive" ! This module's name. if (associated(CS)) then @@ -2113,37 +2110,43 @@ subroutine entrain_diffusive_init(Time, G, GV, US, param_file, diag, CS) endif allocate(CS) + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + CS%diag => diag CS%bulkmixedlayer = (GV%nkml > 0) ! Set default, read and log parameters - call log_version(param_file, mdl, version, "") + if (.not.just_read) call log_version(param_file, mdl, version, "") call get_param(param_file, mdl, "CORRECT_DENSITY", CS%correct_density, & "If true, and USE_EOS is true, the layer densities are "//& "restored toward their target values by the diapycnal "//& "mixing, as described in Hallberg (MWR, 2000).", & - default=.true.) + default=.true., do_not_log=just_read) call get_param(param_file, mdl, "MAX_ENT_IT", CS%max_ent_it, & "The maximum number of iterations that may be used to "//& - "calculate the interior diapycnal entrainment.", default=5) + "calculate the interior diapycnal entrainment.", default=5, do_not_log=just_read) ! In this module, KD is only used to set the default for TOLERANCE_ENT. [m2 s-1] call get_param(param_file, mdl, "KD", Kd, fail_if_missing=.true.) call get_param(param_file, mdl, "DT", dt, & "The (baroclinic) dynamics time step.", units = "s", & - fail_if_missing=.true.) + fail_if_missing=.true., do_not_log=just_read) ! CS%Tolerance_Ent = MAX(100.0*GV%Angstrom_H,1.0e-4*sqrt(dt*Kd)) ! call get_param(param_file, mdl, "TOLERANCE_ENT", CS%Tolerance_Ent, & "The tolerance with which to solve for entrainment values.", & - units="m", default=MAX(100.0*GV%Angstrom_m,1.0e-4*sqrt(dt*Kd)), scale=GV%m_to_H) + units="m", default=MAX(100.0*GV%Angstrom_m,1.0e-4*sqrt(dt*Kd)), scale=GV%m_to_H, do_not_log=just_read) CS%Rho_sig_off = 1000.0*US%kg_m3_to_R - CS%id_Kd = register_diag_field('ocean_model', 'Kd_effective', diag%axesTL, Time, & - 'Diapycnal diffusivity as applied', 'm2 s-1', conversion=US%Z2_T_to_m2_s) - CS%id_diff_work = register_diag_field('ocean_model', 'diff_work', diag%axesTi, Time, & - 'Work actually done by diapycnal diffusion across each interface', & - 'W m-2', conversion=US%RZ3_T3_to_W_m2) + if (.not.just_read) then + CS%id_Kd = register_diag_field('ocean_model', 'Kd_effective', diag%axesTL, Time, & + 'Diapycnal diffusivity as applied', 'm2 s-1', conversion=US%Z2_T_to_m2_s) + CS%id_diff_work = register_diag_field('ocean_model', 'diff_work', diag%axesTi, Time, & + 'Work actually done by diapycnal diffusion across each interface', & + 'W m-2', conversion=US%RZ3_T3_to_W_m2) + endif + + if (just_read) deallocate(CS) end subroutine entrain_diffusive_init diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index 27b316e144..801e573023 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -509,8 +509,7 @@ logical function tidal_mixing_init(Time, G, GV, US, param_file, diag, CS) filename) call safe_alloc_ptr(CS%TKE_Niku,is,ie,js,je) ; CS%TKE_Niku(:,:) = 0.0 call MOM_read_data(filename, 'TKE_input', CS%TKE_Niku, G%domain, timelevel=1, & ! ??? timelevel -aja - scale=US%W_m2_to_RZ3_T3) - CS%TKE_Niku(:,:) = Niku_scale * CS%TKE_Niku(:,:) + scale=Niku_scale*US%W_m2_to_RZ3_T3) call get_param(param_file, mdl, "GAMMA_NIKURASHIN",CS%Gamma_lee, & "The fraction of the lee wave energy that is dissipated "//& @@ -590,7 +589,7 @@ logical function tidal_mixing_init(Time, G, GV, US, param_file, diag, CS) if (CS%use_CVMix_tidal) then CS%id_N2_int = register_diag_field('ocean_model','N2_int',diag%axesTi,Time, & - 'Bouyancy frequency squared, at interfaces', 's-2') !###, conversion=US%s_to_T**2) + 'Bouyancy frequency squared, at interfaces', 's-2', conversion=US%s_to_T**2) !> TODO: add units CS%id_Simmons_coeff = register_diag_field('ocean_model','Simmons_coeff',diag%axesT1,Time, & 'time-invariant portion of the tidal mixing coefficient using the Simmons', '') @@ -781,7 +780,7 @@ subroutine calculate_CVMix_tidal(h, j, G, GV, US, CS, N2_int, Kd_lay, Kv) ! XXX: Temporary de-scaling of N2_int(i,:) into a temporary variable - do k = 1,G%ke+1 + do k=1,G%ke+1 N2_int_i(k) = US%s_to_T**2 * N2_int(i,k) enddo @@ -876,7 +875,7 @@ subroutine calculate_CVMix_tidal(h, j, G, GV, US, CS, N2_int, Kd_lay, Kv) CVmix_tidal_params_user = CS%CVMix_tidal_params) ! XXX: Temporary de-scaling of N2_int(i,:) into a temporary variable - do k = 1,G%ke+1 + do k=1,G%ke+1 N2_int_i(k) = US%s_to_T**2 * N2_int(i,k) enddo diff --git a/src/tracer/ISOMIP_tracer.F90 b/src/tracer/ISOMIP_tracer.F90 index c9bf98f7ff..5503287c50 100644 --- a/src/tracer/ISOMIP_tracer.F90 +++ b/src/tracer/ISOMIP_tracer.F90 @@ -287,7 +287,7 @@ subroutine ISOMIP_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, G if (.not.associated(CS)) return - melt(:,:) = fluxes%iceshelf_melt + melt(:,:) = fluxes%iceshelf_melt(:,:) ! max. melt mmax = MAXVAL(melt(is:ie,js:je)) diff --git a/src/tracer/MOM_lateral_boundary_diffusion.F90 b/src/tracer/MOM_lateral_boundary_diffusion.F90 index 8b9be533d5..f244931376 100644 --- a/src/tracer/MOM_lateral_boundary_diffusion.F90 +++ b/src/tracer/MOM_lateral_boundary_diffusion.F90 @@ -171,16 +171,14 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) ! for diagnostics if (tracer%id_lbdxy_conc > 0 .or. tracer%id_lbdxy_cont > 0 .or. tracer%id_lbdxy_cont_2d > 0) then - tendency(:,:,:) = 0.0 + tendency(:,:,:) = 0.0 endif - do j = G%jsc-1, G%jec+1 - ! Interpolate state to interface - do i = G%isc-1, G%iec+1 - call build_reconstructions_1d( CS%remap_CS, G%ke, h(i,j,:), tracer%t(i,j,:), ppoly0_coefs(i,j,:,:), & - ppoly0_E(i,j,:,:), ppoly_S, remap_method, GV%H_subroundoff, GV%H_subroundoff) - enddo - enddo + ! Interpolate state to interface + do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1 + call build_reconstructions_1d( CS%remap_CS, G%ke, h(i,j,:), tracer%t(i,j,:), ppoly0_coefs(i,j,:,:), & + ppoly0_E(i,j,:,:), ppoly_S, remap_method, GV%H_subroundoff, GV%H_subroundoff) + enddo ; enddo ! Diffusive fluxes in the i-direction uFlx(:,:,:) = 0. vFlx(:,:,:) = 0. @@ -253,41 +251,41 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) if (tracer%id_lbd_dfy>0) call post_data(tracer%id_lbd_dfy, vFlx*Idt, CS%diag) if (tracer%id_lbd_dfx_2d>0) then uwork_2d(:,:) = 0. - do k=1,GV%ke; do j=G%jsc,G%jec; do I=G%isc-1,G%iec + do k=1,GV%ke ; do j=G%jsc,G%jec ; do I=G%isc-1,G%iec uwork_2d(I,j) = uwork_2d(I,j) + (uFlx(I,j,k) * Idt) - enddo; enddo; enddo + enddo ; enddo ; enddo call post_data(tracer%id_lbd_dfx_2d, uwork_2d, CS%diag) endif if (tracer%id_lbd_dfy_2d>0) then vwork_2d(:,:) = 0. - do k=1,GV%ke; do J=G%jsc-1,G%jec; do i=G%isc,G%iec + do k=1,GV%ke ; do J=G%jsc-1,G%jec ; do i=G%isc,G%iec vwork_2d(i,J) = vwork_2d(i,J) + (vFlx(i,J,k) * Idt) - enddo; enddo; enddo + enddo ; enddo ; enddo call post_data(tracer%id_lbd_dfy_2d, vwork_2d, CS%diag) endif ! post tendency of tracer content if (tracer%id_lbdxy_cont > 0) then - call post_data(tracer%id_lbdxy_cont, tendency(:,:,:), CS%diag) + call post_data(tracer%id_lbdxy_cont, tendency, CS%diag) endif ! post depth summed tendency for tracer content if (tracer%id_lbdxy_cont_2d > 0) then tendency_2d(:,:) = 0. - do j = G%jsc,G%jec ; do i = G%isc,G%iec - do k = 1, GV%ke + do j=G%jsc,G%jec ; do i=G%isc,G%iec + do k=1,GV%ke tendency_2d(i,j) = tendency_2d(i,j) + tendency(i,j,k) enddo enddo ; enddo - call post_data(tracer%id_lbdxy_cont_2d, tendency_2d(:,:), CS%diag) + call post_data(tracer%id_lbdxy_cont_2d, tendency_2d, CS%diag) endif ! post tendency of tracer concentration; this step must be ! done after posting tracer content tendency, since we alter ! the tendency array and its units. if (tracer%id_lbdxy_conc > 0) then - do k = 1, GV%ke ; do j = G%jsc,G%jec ; do i = G%isc,G%iec + do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec tendency(i,j,k) = tendency(i,j,k) / ( h(i,j,k) + GV%H_subroundoff ) enddo ; enddo ; enddo call post_data(tracer%id_lbdxy_conc, tendency, CS%diag) @@ -306,9 +304,9 @@ real function bulk_average(boundary, nk, deg, h, hBLT, phi, ppoly0_E, ppoly0_coe real, dimension(nk) :: h !< Layer thicknesses [H ~> m or kg m-2] real :: hBLT !< Depth of the boundary layer [H ~> m or kg m-2] real, dimension(nk) :: phi !< Scalar quantity - real, dimension(nk,2) :: ppoly0_E(:,:) !< Edge value of polynomial - real, dimension(nk,deg+1) :: ppoly0_coefs(:,:) !< Coefficients of polynomial - integer :: method !< Remapping scheme to use + real, dimension(nk,2) :: ppoly0_E !< Edge value of polynomial + real, dimension(nk,deg+1) :: ppoly0_coefs !< Coefficients of polynomial + integer :: method !< Remapping scheme to use integer :: k_top !< Index of the first layer within the boundary real :: zeta_top !< Fraction of the layer encompassed by the bottom boundary layer diff --git a/src/tracer/MOM_offline_main.F90 b/src/tracer/MOM_offline_main.F90 index 65f83ecfea..b7af9849b3 100644 --- a/src/tracer/MOM_offline_main.F90 +++ b/src/tracer/MOM_offline_main.F90 @@ -28,7 +28,7 @@ module MOM_offline_main use MOM_offline_aux, only : distribute_residual_uh_upwards, distribute_residual_vh_upwards use MOM_opacity, only : opacity_CS, optics_type use MOM_open_boundary, only : ocean_OBC_type -use MOM_time_manager, only : time_type +use MOM_time_manager, only : time_type, real_to_time use MOM_tracer_advect, only : tracer_advect_CS, advect_tracer use MOM_tracer_diabatic, only : applyTracerBoundaryFluxesInOut use MOM_tracer_flow_control, only : tracer_flow_control_CS, call_tracer_column_fns, call_tracer_stocks @@ -79,7 +79,8 @@ module MOM_offline_main integer :: start_index !< Timelevel to start integer :: iter_no !< Timelevel to start integer :: numtime !< How many timelevels in the input fields - integer :: accumulated_time !< Length of time accumulated in the current offline interval + type(time_type) :: accumulated_time !< Length of time accumulated in the current offline interval + type(time_type) :: vertical_time !< The next value of accumulate_time at which to apply vertical processes ! Index of each of the variables to be read in with separate indices for each variable if they ! are set off from each other in time integer :: ridx_sum = -1 !< Read index offset of the summed variables @@ -722,9 +723,9 @@ subroutine offline_diabatic_ale(fluxes, Time_start, Time_end, CS, h_pre, eatr, e ! Add diurnal cycle for shortwave radiation (only used if run in ocean-only mode) if (CS%diurnal_SW .and. CS%read_sw) then - sw(:,:) = fluxes%sw - sw_vis(:,:) = fluxes%sw_vis_dir - sw_nir(:,:) = fluxes%sw_nir_dir + sw(:,:) = fluxes%sw(:,:) + sw_vis(:,:) = fluxes%sw_vis_dir(:,:) + sw_nir(:,:) = fluxes%sw_nir_dir(:,:) call offline_add_diurnal_SW(fluxes, CS%G, Time_start, Time_end) endif @@ -738,9 +739,9 @@ subroutine offline_diabatic_ale(fluxes, Time_start, Time_end, CS, h_pre, eatr, e CS%G, CS%GV, CS%US, CS%tv, CS%optics, CS%tracer_flow_CSp, CS%debug) if (CS%diurnal_SW .and. CS%read_sw) then - fluxes%sw(:,:) = sw - fluxes%sw_vis_dir(:,:) = sw_vis - fluxes%sw_nir_dir(:,:) = sw_nir + fluxes%sw(:,:) = sw(:,:) + fluxes%sw_vis_dir(:,:) = sw_vis(:,:) + fluxes%sw_nir_dir(:,:) = sw_nir(:,:) endif if (CS%debug) then @@ -1200,7 +1201,7 @@ end subroutine post_offline_convergence_diags !> Extracts members of the offline main control structure. All arguments are optional except !! the control structure itself -subroutine extract_offline_main(CS, uhtr, vhtr, eatr, ebtr, h_end, accumulated_time, & +subroutine extract_offline_main(CS, uhtr, vhtr, eatr, ebtr, h_end, accumulated_time, vertical_time, & dt_offline, dt_offline_vertical, skip_diffusion) type(offline_transport_CS), target, intent(in ) :: CS !< Offline control structure ! Returned optional arguments @@ -1212,9 +1213,10 @@ subroutine extract_offline_main(CS, uhtr, vhtr, eatr, ebtr, h_end, accumulated_t !! one time step [H ~> m or kg m-2] real, dimension(:,:,:), optional, pointer :: h_end !< Thicknesses at the end of offline timestep !! [H ~> m or kg m-2] - !### Why are the following variables integers? - integer, optional, pointer :: accumulated_time !< Length of time accumulated in the - !! current offline interval [s] + type(time_type), optional, pointer :: accumulated_time !< Length of time accumulated in the + !! current offline interval + type(time_type), optional, pointer :: vertical_time !< The next value of accumulate_time at which to + !! vertical processes real, optional, intent( out) :: dt_offline !< Timestep used for offline tracers [T ~> s] real, optional, intent( out) :: dt_offline_vertical !< Timestep used for calls to tracer !! vertical physics [T ~> s] @@ -1229,6 +1231,7 @@ subroutine extract_offline_main(CS, uhtr, vhtr, eatr, ebtr, h_end, accumulated_t ! Pointers to integer members which need to be modified if (present(accumulated_time)) accumulated_time => CS%accumulated_time + if (present(vertical_time)) vertical_time => CS%vertical_time ! Return value of non-modified integers if (present(dt_offline)) dt_offline = CS%dt_offline @@ -1414,7 +1417,8 @@ subroutine offline_transport_init(param_file, CS, diabatic_CSp, G, GV, US) end select ! Set the accumulated time to zero - CS%accumulated_time = 0 + CS%accumulated_time = real_to_time(0.0) + CS%vertical_time = CS%accumulated_time ! Set the starting read index for time-averaged and snapshotted fields CS%ridx_sum = CS%start_index if (CS%fields_are_offset) CS%ridx_snap = next_modulo_time(CS%start_index,CS%numtime)