diff --git a/.circleci/config.yml b/.circleci/config.yml old mode 100644 new mode 100755 diff --git a/.gitignore b/.gitignore old mode 100644 new mode 100755 diff --git a/.gitmodules b/.gitmodules old mode 100644 new mode 100755 diff --git a/cime_config/machines/cmake_macros/gnu_chicoma-cpu.cmake b/cime_config/machines/cmake_macros/gnu_chicoma-cpu.cmake index b6235f57ab26..0fcef4651937 100644 --- a/cime_config/machines/cmake_macros/gnu_chicoma-cpu.cmake +++ b/cime_config/machines/cmake_macros/gnu_chicoma-cpu.cmake @@ -23,3 +23,5 @@ set(MPIFC "ftn") set(SCC "gcc") set(SCXX "g++") set(SFC "gfortran") +set(PETSC_PATH "$ENV{PETSC_PATH}") +string(APPEND SLIBS " -L$ENV{PETSC_PATH}/lib -lpetsc") diff --git a/cime_config/machines/cmake_macros/intel_anvil.cmake b/cime_config/machines/cmake_macros/intel_anvil.cmake index 290869f36c80..5745f28563e3 100644 --- a/cime_config/machines/cmake_macros/intel_anvil.cmake +++ b/cime_config/machines/cmake_macros/intel_anvil.cmake @@ -39,3 +39,5 @@ endif() set(NETCDF_C_PATH "$ENV{NETCDF_C_PATH}") set(NETCDF_FORTRAN_PATH "$ENV{NETCDF_FORTRAN_PATH}") set(PNETCDF_PATH "$ENV{PNETCDF_PATH}") +set(PETSC_PATH "$ENV{PETSC_PATH}") +string(APPEND SLIBS " -L$ENV{PETSC_PATH}/lib -lpetsc") diff --git a/cime_config/machines/config_machines.xml b/cime_config/machines/config_machines.xml index b1e899391d21..1ae9b0390254 100644 --- a/cime_config/machines/config_machines.xml +++ b/cime_config/machines/config_machines.xml @@ -2335,6 +2335,7 @@ netcdf-cxx/4.2-gkqc6fq netcdf-fortran/4.4.4-eanrh5t parallel-netcdf/1.11.0-y3nmmej + petsc-3.16.1-intel-20.0.4-3vcputc openmpi/4.1.1-v3b3npd @@ -2377,6 +2378,7 @@ $SHELL{dirname $(dirname $(which nc-config))} $SHELL{dirname $(dirname $(which nf-config))} $SHELL{dirname $(dirname $(which pnetcdf_version))} + /lcrc/soft/climate/compass/anvil/spack/spack_for_mache_1.4.1/var/spack/environments/compass_1_1_0_intel_impi_netlib_lapack_petsc/.spack-env/view /lcrc/group/e3sm/soft/perl/5.26.0/bin:$ENV{PATH} @@ -4065,6 +4067,10 @@ software MPI_Bcast + + /usr/projects/climate/SHARED_CLIMATE/software/chicoma-cpu/petsc/gcc-12.1.0/cray-mpich-8.1.21 + /usr/projects/climate/SHARED_CLIMATE/software/chicoma-cpu/petsc/gcc-12.1.0/cray-mpich-8.1.21/lib:$ENV{LD_LIBRARY_PATH} + -1 diff --git a/components/mpas-framework/Makefile b/components/mpas-framework/Makefile index 06865333cdd3..7795b9ac9c36 100644 --- a/components/mpas-framework/Makefile +++ b/components/mpas-framework/Makefile @@ -508,7 +508,26 @@ CPPINCLUDES = FCINCLUDES = LIBS = -# +# Add links to PETSC links if requested +ifeq "$(USE_PETSC)" "true" +ifndef PETSC +$(error PETSC is not set. Please set PETSC to the PETSC install directory when USE_PETSC=true) +endif +ifneq (, $(shell ls $(PETSC)/lib/libpetsc.*)) + LIBS += -L$(PETSC)/lib + CPPINCLUDES += -I$(PETSC)/include + FCINCLUDES += -I$(PETSC)/include +else ifneq (, $(shell ls $(PETSC)/libpetsc.*)) + LIBS += -L$(PETSC) + CPPINCLUDES += -I$(PETSC) + FCINCLUDES += -I$(PETSC) +else +$(error libpetsc.* does NOT exist in $(PETSC) or $(PETSC)/lib) +endif + LIBS += -lpetsc + override CPPFLAGS += -DUSE_PETSC +endif + # If user has indicated a PIO2 library, define USE_PIO2 pre-processor macro # ifeq "$(USE_PIO2)" "true" @@ -524,7 +543,7 @@ ifneq ($(wildcard $(PIO)/lib), ) else PIO_LIB = $(PIO) endif -LIBS = -L$(PIO_LIB) +LIBS += -L$(PIO_LIB) # # Regardless of PIO library version, look for an include subdirectory of PIO path diff --git a/components/mpas-ocean/bld/build-namelist b/components/mpas-ocean/bld/build-namelist index 90967685be19..3603af9b337f 100755 --- a/components/mpas-ocean/bld/build-namelist +++ b/components/mpas-ocean/bld/build-namelist @@ -526,6 +526,8 @@ add_default($nl, 'config_use_mom_del2'); add_default($nl, 'config_mom_del2'); add_default($nl, 'config_use_tracer_del2'); add_default($nl, 'config_tracer_del2'); +add_default($nl, 'config_use_vertMom_del2'); +add_default($nl, 'config_vertMom_del2'); ############################# # Namelist group: hmix_del4 # @@ -836,6 +838,19 @@ add_default($nl, 'config_vegetation_drag_coefficient'); add_default($nl, 'config_density0'); +########################################### +# Namelist group: nonhydrostatic_pressure # +########################################### + +add_default($nl, 'config_enable_nonhydrostatic_mode'); +add_default($nl, 'config_nonhydrostatic_solve_surface_boundary_condition'); +add_default($nl, 'config_nonhydrostatic_remove_rhs_mean'); +add_default($nl, 'config_nonhydrostatic_preconditioner'); +add_default($nl, 'config_nonhydrostatic_solver_type'); +add_default($nl, 'config_petsc_rtol'); +add_default($nl, 'config_petsc_atol'); +add_default($nl, 'config_petsc_maxit'); + ##################################### # Namelist group: pressure_gradient # ##################################### @@ -956,6 +971,9 @@ add_default($nl, 'config_disable_vel_topographic_wave_drag'); add_default($nl, 'config_disable_vel_explicit_bottom_drag'); add_default($nl, 'config_disable_vel_vmix'); add_default($nl, 'config_disable_vel_vadv'); +add_default($nl, 'config_disable_vvel_all_tend'); +add_default($nl, 'config_disable_vvel_coriolis'); +add_default($nl, 'config_disable_vvel_adv'); add_default($nl, 'config_disable_tr_all_tend'); add_default($nl, 'config_disable_tr_adv'); add_default($nl, 'config_disable_tr_hmix'); @@ -1744,6 +1762,7 @@ my @groups = qw(run_modes rayleigh_damping vegetation_drag ocean_constants + nonhydrostatic_pressure pressure_gradient eos eos_linear diff --git a/components/mpas-ocean/bld/build-namelist-group-list b/components/mpas-ocean/bld/build-namelist-group-list index 7d5b8cdd473f..7529f6dec309 100644 --- a/components/mpas-ocean/bld/build-namelist-group-list +++ b/components/mpas-ocean/bld/build-namelist-group-list @@ -26,6 +26,7 @@ my @groups = qw(run_modes rayleigh_damping vegetation_drag ocean_constants + nonhydrostatic_pressure pressure_gradient eos eos_linear diff --git a/components/mpas-ocean/bld/build-namelist-section b/components/mpas-ocean/bld/build-namelist-section index 8dcdab06d6e4..1701fbbe3815 100644 --- a/components/mpas-ocean/bld/build-namelist-section +++ b/components/mpas-ocean/bld/build-namelist-section @@ -65,6 +65,8 @@ add_default($nl, 'config_use_mom_del2'); add_default($nl, 'config_mom_del2'); add_default($nl, 'config_use_tracer_del2'); add_default($nl, 'config_tracer_del2'); +add_default($nl, 'config_use_vertMom_del2'); +add_default($nl, 'config_vertMom_del2'); ############################# # Namelist group: hmix_del4 # @@ -355,6 +357,19 @@ add_default($nl, 'config_vegetation_drag_coefficient'); add_default($nl, 'config_density0'); +########################################### +# Namelist group: nonhydrostatic_pressure # +########################################### + +add_default($nl, 'config_enable_nonhydrostatic_mode'); +add_default($nl, 'config_nonhydrostatic_solve_surface_boundary_condition'); +add_default($nl, 'config_nonhydrostatic_remove_rhs_mean'); +add_default($nl, 'config_nonhydrostatic_preconditioner'); +add_default($nl, 'config_nonhydrostatic_solver_type'); +add_default($nl, 'config_petsc_rtol'); +add_default($nl, 'config_petsc_atol'); +add_default($nl, 'config_petsc_maxit'); + ##################################### # Namelist group: pressure_gradient # ##################################### @@ -475,6 +490,9 @@ add_default($nl, 'config_disable_vel_topographic_wave_drag'); add_default($nl, 'config_disable_vel_explicit_bottom_drag'); add_default($nl, 'config_disable_vel_vmix'); add_default($nl, 'config_disable_vel_vadv'); +add_default($nl, 'config_disable_vvel_all_tend'); +add_default($nl, 'config_disable_vvel_coriolis'); +add_default($nl, 'config_disable_vvel_adv'); add_default($nl, 'config_disable_tr_all_tend'); add_default($nl, 'config_disable_tr_adv'); add_default($nl, 'config_disable_tr_hmix'); diff --git a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml index d2cbcaee3630..bbad055eeb66 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml +++ b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml @@ -96,6 +96,8 @@ 1000.0 .false. 10.0 +.false. +10.0 .true. @@ -395,6 +397,16 @@ 1026.0 + +.false. +'noGradient' +.true. +'jacobi' +'pipecgrr' +1.0e-10 +1.0E-20 +10000 + 'Jacobian_from_TS' 0.5 @@ -509,6 +521,9 @@ .false. .false. .false. +.false. +.false. +.false. .false. .false. .false. diff --git a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml index 2b1b736252aa..953ab0fb15e4 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml +++ b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml @@ -286,6 +286,22 @@ Valid values: any positive real Default: Defined in namelist_defaults.xml + +If true, Laplacian horizontal mixing is used on the vertical momentum equation. + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +Horizontal diffusion. + +Valid values: any positive real +Default: Defined in namelist_defaults.xml + + @@ -1821,6 +1837,71 @@ Valid values: any positive real, but typically 1000-1035 Default: Defined in namelist_defaults.xml + + + +Flag to enable the nonhydrostatic pressure calculation + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +setting for surface boundary condition for the pressure correction solve + +Valid values: noGradient or no pressure +Default: Defined in namelist_defaults.xml + + + +flag to enabling mean from the right hand side of the nonhydro solve + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +String describing preconditioner to use, recommended options -- jacobi, bjacobi, ilu + +Valid values: jacobi, bjacobi, ilu +Default: Defined in namelist_defaults.xml + + + +String describing the solver context, recommended types -- pipecgrr, pipecg, cg, bcgs, gmres + +Valid values: pipecgrr, pipecg, cg, bcgs, gmres +Default: Defined in namelist_defaults.xml + + + +relative tolerance for the norm in petsc convergence + +Valid values: very small positive reals +Default: Defined in namelist_defaults.xml + + + +absolute tolerance for the norm in petsc convergence + +Valid values: very small positive reals +Default: Defined in namelist_defaults.xml + + + +maximum number of iterations in petsc solve + +Valid values: positive integers +Default: Defined in namelist_defaults.xml + @@ -2419,6 +2500,30 @@ Valid values: .true. or .false. Default: Defined in namelist_defaults.xml + +Disables all tendencies on vertical momentum field. + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +Disables tendencies on the vetical velocity field from the Coriolis force. + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +Disables tendencies on the vertical velocity field from vertical and horizontal advection. + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + Disables all tendencies on tracer fields. diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index 467d831ff299..a3513d2f1080 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -237,6 +237,14 @@ description="If true, Laplacian horizontal mixing is used on the tracer equation." possible_values=".true. or .false." /> + + + + + + + + + + + + + + + + + + @@ -1877,14 +1940,16 @@ - - - - - - - - + + + + + + + + + + + + + + diff --git a/components/mpas-ocean/src/driver/mpas_ocn_core_interface.F b/components/mpas-ocean/src/driver/mpas_ocn_core_interface.F index 3359d68d516f..ff8d60280e35 100644 --- a/components/mpas-ocean/src/driver/mpas_ocn_core_interface.F +++ b/components/mpas-ocean/src/driver/mpas_ocn_core_interface.F @@ -110,6 +110,7 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{ integer :: err_tmp logical, pointer :: forwardModeActive + logical, pointer :: nonhydrostaticPKGActive logical, pointer :: analysisModeActive logical, pointer :: initModeActive logical, pointer :: thicknessFilterActive @@ -165,6 +166,7 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{ logical, pointer :: config_use_GM logical, pointer :: config_submesoscale_enable logical, pointer :: config_use_Redi + logical, pointer :: config_enable_nonhydrostatic_mode logical, pointer :: config_use_vegetation_drag logical, pointer :: config_use_time_varying_atmospheric_forcing logical, pointer :: config_use_time_varying_land_ice_forcing @@ -349,6 +351,15 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{ topographicWaveDragPKGActive = .true. end if + ! + ! test for nonHydrostatic pressure + ! + call mpas_pool_get_package(packagePool, 'nonhydrostaticPKGActive', nonhydrostaticPKGActive) + call mpas_pool_get_config(configPool, 'config_enable_nonhydrostatic_mode', config_enable_nonhydrostatic_mode) + if (config_enable_nonhydrostatic_mode) then + nonhydrostaticPKGActive = .true. + end if + ! ! test for use of gm ! diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F index 0e5235791f14..66d196a92fbc 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F @@ -49,6 +49,12 @@ module ocn_forward_mode use ocn_thick_ale use ocn_thick_surface_flux + use ocn_vvel_pgrad + use ocn_vvel_coriolis + use ocn_vvel_advection + use ocn_vvel_hmix_del2 + use ocn_vvel_hmix_del4 + use ocn_vel_pressure_grad use ocn_vel_vadv use ocn_vel_hmix @@ -76,6 +82,7 @@ module ocn_forward_mode use ocn_tracer_CFC use ocn_tracer_surface_restoring use ocn_gm + use ocn_compute_nonhydrostatic_pressure use ocn_submesoscale_eddies use ocn_stokes_drift @@ -310,6 +317,21 @@ function ocn_forward_mode_init(domain, startTimeStamp) result(ierr)!{{{ timeStep = mpas_get_clock_timestep(domain % clock, ierr=err_tmp) call mpas_get_timeInterval(timeStep, dt=dt) + if(config_enable_nonhydrostatic_mode) then + call ocn_vvel_advection_init(err_tmp) + ierr = ior(ierr, err_tmp) + call ocn_vvel_coriolis_init(err_tmp) + ierr = ior(ierr, err_tmp) + call ocn_vvel_pgf_init(err_tmp) + ierr = ior(ierr, err_tmp) + call ocn_vvel_hmix_del2_init(err_tmp) + ierr = ior(ierr, err_tmp) + call ocn_vvel_hmix_del4_init(err_tmp) + ierr = ior(ierr, err_tmp) + call ocn_nonhydrostatic_solver_init(domain % dminfo % comm, err_tmp) + ierr = ior(ierr, err_tmp) + endif + call ocn_init_routines_block(block, dt, ierr) if(ierr.eq.1) then @@ -867,6 +889,10 @@ function ocn_forward_mode_finalize(domain) result(iErr)!{{{ call mpas_decomp_destroy_decomp_list(domain % decompositions) + if(config_enable_nonhydrostatic_mode) then + call ocn_nonhydrostatic_solver_finalize() + end if + end function ocn_forward_mode_finalize!}}} end module ocn_forward_mode diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F index 58755c871d12..ff5aaa0d535c 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F @@ -35,6 +35,7 @@ module ocn_time_integration_rk4 use ocn_gm use ocn_submesoscale_eddies + use ocn_compute_nonhydrostatic_pressure use ocn_equation_of_state use ocn_vmix use ocn_time_average_coupled @@ -139,6 +140,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ logical, pointer :: config_verify_not_dry logical, pointer :: config_prevent_drying logical, pointer :: config_zero_drying_velocity + logical, pointer :: config_enable_nonhydrostatic_mode character (len=StrKIND), pointer :: config_land_ice_flux_mode ! Diagnostics Indices @@ -151,18 +153,22 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ real (kind=RKIND), dimension(:,:), pointer :: normalVelocityProvis, layerThicknessProvis real (kind=RKIND), dimension(:,:), pointer :: highFreqThicknessProvis real (kind=RKIND), dimension(:,:), pointer :: lowFreqDivergenceProvis + real (kind=RKIND), dimension(:,:), pointer :: vertVelocityNonhydroProvis real (kind=RKIND), dimension(:,:,:), pointer :: tracersGroupProvis ! Tend Array Pointers real (kind=RKIND), dimension(:,:), pointer :: highFreqThicknessTend, lowFreqDivergenceTend, normalVelocityTend, & - layerThicknessTend + layerThicknessTend, vertVelocityNonhydroTend real (kind=RKIND), dimension(:,:,:), pointer :: tracersGroupTend ! State Array Pointers real (kind=RKIND), dimension(:,:), pointer :: normalVelocityCur, normalVelocityNew + real (kind=RKIND), dimension(:,:), pointer :: vertVelocityNonhydroCur, vertVelocityNonhydroNew real (kind=RKIND), dimension(:,:), pointer :: layerThicknessCur, layerThicknessNew real (kind=RKIND), dimension(:,:), pointer :: highFreqThicknessCur, highFreqThicknessNew real (kind=RKIND), dimension(:,:), pointer :: lowFreqDivergenceCur, lowFreqDivergenceNew + real (kind=RKIND), dimension(:,:), pointer :: nonhydrostaticPressureCur, nonhydrostaticPressureNew + real (kind=RKIND), dimension(:), pointer :: sshCur, sshNew real (kind=RKIND), dimension(:,:,:), pointer :: tracerGroup, tracersCur, tracersNew @@ -214,6 +220,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ call mpas_pool_get_config(domain % configs, 'config_zero_drying_velocity', config_zero_drying_velocity) call mpas_pool_get_config(domain % configs, 'config_use_tidal_forcing', config_use_tidal_forcing) call mpas_pool_get_config(domain % configs, 'config_tidal_forcing_type', config_tidal_forcing_type) + call mpas_pool_get_config(domain % configs, 'config_enable_nonhydrostatic_mode', config_enable_nonhydrostatic_mode) ! ! Initialize time_levs(2) with state at current time @@ -241,6 +248,11 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ call mpas_pool_get_array(statePool, 'layerThickness', layerThicknessCur, 1) call mpas_pool_get_array(statePool, 'layerThickness', layerThicknessNew, 2) + if(config_enable_nonhydrostatic_mode) then + call mpas_pool_get_array(statePool, 'vertVelocityNonhydro', vertVelocityNonhydroCur, 1) + call mpas_pool_get_array(statePool, 'vertVelocityNonhydro', vertVelocityNonhydroNew, 2) + end if + call mpas_pool_get_array(statePool, 'highFreqThickness', highFreqThicknessCur, 1) call mpas_pool_get_array(statePool, 'highFreqThickness', highFreqThicknessNew, 2) call mpas_pool_get_array(statePool, 'lowFreqDivergence', lowFreqDivergenceCur, 1) @@ -270,6 +282,18 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ !$omp end do !$omp end parallel + if(config_enable_nonhydrostatic_mode) then + !$omp parallel + !$omp do schedule(runtime) private(k) + do iCell = 1, nCells + do k = 1, maxLevelCell(iCell) + vertVelocityNonhydroNew(k,iCell) = vertVelocityNonhydroCur(k, iCell) + end do + end do + !$omp end do + !$omp end parallel + end if + call mpas_pool_begin_iteration(tracersPool) do while ( mpas_pool_get_next_member(tracersPool, groupItr) ) @@ -508,12 +532,22 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ call mpas_timer_stop("RK4 vel/thick tendency computations") - ! Update halos for prognostic variables. + if (config_enable_nonhydrostatic_mode) then + call mpas_timer_start("RK4 vvel tendency computations") + block => domain % blocklist + do while (associated(block)) + call ocn_time_integrator_rk4_compute_vvel_tends( block, dt, rk_substep_weights(rk_step), err ) + block => block % next + end do + call mpas_timer_stop("RK4 vvel tendency computations") + end if + ! Update halos for prognostic variables. call mpas_timer_start("RK4 vel/thick prognostic halo update") call mpas_dmpar_field_halo_exch(domain, 'tendNormalVelocity') call mpas_dmpar_field_halo_exch(domain, 'tendLayerThickness') + call mpas_dmpar_field_halo_exch(domain, 'tendVertVelocityNonhydro') call mpas_timer_stop("RK4 vel/thick prognostic halo update") @@ -603,7 +637,8 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ ! Rescale tracers block => domain % blocklist do while(associated(block)) - call ocn_time_integrator_rk4_cleanup(block, dt, err) + dminfo = domain % dminfo + call ocn_time_integrator_rk4_cleanup(block, dt, dminfo, err) block => block % next end do @@ -619,6 +654,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ call mpas_pool_get_subpool(statePool, 'tracers', tracersPool) call mpas_dmpar_field_halo_exch(domain, 'normalVelocity', timeLevel=2) + call mpas_dmpar_field_halo_exch(domain, 'vertVelocityNonhydro', timeLevel=2) call mpas_pool_begin_iteration(tracersPool) do while ( mpas_pool_get_next_member(tracersPool, groupItr) ) @@ -659,6 +695,22 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ call mpas_pool_get_array(tracersPool, 'activeTracers',activeTracersCur, 1) call mpas_pool_get_array(tracersPool, 'activeTracers',activeTracersNew, 2) + if (config_enable_nonhydrostatic_mode) then + + call mpas_pool_get_array(statePool, 'nonhydrostaticPressure', nonhydrostaticPressureNew, 2) + call mpas_pool_get_array(statePool, 'nonhydrostaticPressure', nonhydrostaticPressureCur, 1) + + call mpas_timer_start('RK4 NH pressure tend halo pressure') + call mpas_dmpar_field_halo_exch(domain, 'nonhydrostaticPressure', timeLevel=2) + call mpas_dmpar_field_halo_exch(domain, 'nonhydrostaticPressure', timeLevel=1) + call mpas_dmpar_field_halo_exch(domain, 'nhPressureCorrection') + call mpas_timer_stop('RK4 NH pressure tend halo pressure') + + call ocn_nonhydrostatic_pressure_update_velocity(normalVelocityNew, vertVelocityNonhydroNew, & + layerThicknessNew, layerThickEdgeMean, nonhydrostaticPressureNew, nonhydrostaticPressureCur, & + nhPressureCorrection, dt) + end if + if (config_prescribe_velocity) then !$omp parallel !$omp do schedule(runtime) @@ -667,6 +719,16 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ end do !$omp end do !$omp end parallel + + if (config_enable_nonhydrostatic_mode) then + !$omp parallel + !$omp do schedule(runtime) + do iCell = 1, nCells + vertVelocityNonhydroNew(:, iCell) = vertVelocityNonhydroCur(:, iCell) + end do + !$omp end do + !$omp end parallel + end if end if if (config_disable_tr_all_tend) then @@ -882,6 +944,7 @@ subroutine ocn_time_integrator_rk4_compute_vel_tends(domain, block, dt, & real (kind=RKIND), dimension(:), pointer :: sshCur real (kind=RKIND), dimension(:, :), pointer :: layerThicknessCur, normalVelocityCur + real (kind=RKIND), dimension(:, :), pointer :: vertVelocityNonhydro real (kind=RKIND), dimension(:, :), pointer :: normalVelocityProvis, highFreqThicknessProvis logical, pointer :: config_filter_btr_mode @@ -903,6 +966,7 @@ subroutine ocn_time_integrator_rk4_compute_vel_tends(domain, block, dt, & call mpas_pool_get_array(statePool, 'layerThickness', layerThicknessCur, 1) call mpas_pool_get_array(statePool, 'ssh', sshCur, 1) call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocityCur, 1) + call mpas_pool_get_array(statePool, 'vertVelocityNonhydro', vertVelocityNonhydro, 1) call mpas_pool_get_array(provisStatePool, 'normalVelocity', normalVelocityProvis, 1) call mpas_pool_get_array(provisStatePool, 'highFreqThickness', highFreqThicknessProvis, 1) @@ -911,19 +975,81 @@ subroutine ocn_time_integrator_rk4_compute_vel_tends(domain, block, dt, & if (associated(highFreqThicknessProvis)) then call ocn_vert_transport_velocity_top(verticalMeshPool, & layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, & - sshCur, rkSubstepWeight, & - vertAleTransportTop, err, highFreqThicknessProvis) + sshCur, rkSubstepWeight, vertAleTransportTop, err, & + highFreqThicknessProvis) else call ocn_vert_transport_velocity_top(verticalMeshPool, & layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, & - sshCur, rkSubstepWeight, & - vertAleTransportTop, err) + sshCur, rkSubstepWeight, vertAleTransportTop, err) endif call ocn_tend_vel(domain, tendPool, provisStatePool, forcingPool, 1, dminfo, dt) end subroutine ocn_time_integrator_rk4_compute_vel_tends!}}} + subroutine ocn_time_integrator_rk4_compute_vvel_tends(block, dt, & + rkSubstepWeight, err) + + type (block_type), intent(in) :: block + real (kind=RKIND), intent(in) :: dt + real (kind=RKIND), intent(in) :: rkSubstepWeight + integer, intent(out) :: err + + type (mpas_pool_type), pointer :: meshPool, verticalMeshPool + type (mpas_pool_type), pointer :: statePool, forcingPool + type (mpas_pool_type), pointer :: scratchPool, tendPool, provisStatePool + type (mpas_pool_type), pointer :: tracersPool + + real (kind=RKIND), dimension(:), pointer :: sshCur + real (kind=RKIND), dimension(:, :), pointer :: layerThicknessCur, normalVelocityCur + real (kind=RKIND), dimension(:, :), pointer :: vertVelocityNonhydroCur + real (kind=RKIND), dimension(:, :), pointer :: normalVelocityProvis, highFreqThicknessProvis + real (kind=RKIND), dimension(:, :), pointer :: vertVelocityNonhydroProvis + + integer :: iEdge + integer, pointer :: nEdges + logical, pointer :: config_filter_btr_mode + + err = 0 + + call mpas_pool_get_config(block % configs, 'config_filter_btr_mode', config_filter_btr_mode) + + call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) + call mpas_pool_get_subpool(block % structs, 'verticalMesh', verticalMeshPool) + call mpas_pool_get_subpool(block % structs, 'state', statePool) + call mpas_pool_get_subpool(block % structs, 'forcing', forcingPool) + call mpas_pool_get_subpool(block % structs, 'scratch', scratchPool) + call mpas_pool_get_subpool(block % structs, 'tend', tendPool) + call mpas_pool_get_subpool(block % structs, 'provis_state', provisStatePool) + + call mpas_pool_get_subpool(statePool, 'tracers', tracersPool) + + call mpas_pool_get_dimension(block % dimensions, 'nEdges', nEdges) + call mpas_pool_get_array(statePool, 'layerThickness', layerThicknessCur, 1) + call mpas_pool_get_array(statePool, 'ssh', sshCur, 1) + call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocityCur, 1) + call mpas_pool_get_array(statePool, 'vertVelocityNonhydro', vertVelocityNonhydroCur, 1) + + call mpas_pool_get_array(provisStatePool, 'vertVelocityNonhydro', vertVelocityNonhydroProvis, 1) + call mpas_pool_get_array(provisStatePool, 'normalVelocity', normalVelocityProvis, 1) + call mpas_pool_get_array(provisStatePool, 'highFreqThickness', highFreqThicknessProvis, 1) + + ! advection of u uses u, while advection of layerThickness and tracers use normalTransportVelocity. + if (associated(highFreqThicknessProvis)) then + call ocn_vert_transport_velocity_top(verticalMeshPool, & + layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, & + sshCur, rkSubstepWeight, vertAleTransportTop, err, & + highFreqThicknessProvis) + else + call ocn_vert_transport_velocity_top(verticalMeshPool, & + layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, & + sshCur, rkSubstepWeight, vertAleTransportTop, err) + endif + + call ocn_tend_vvel(tendPool, provisStatePool, forcingPool, meshPool, 1, dt) + + end subroutine ocn_time_integrator_rk4_compute_vvel_tends + subroutine ocn_time_integrator_rk4_compute_thick_tends(block, dt, rkSubstepWeight, err)!{{{ type (block_type), intent(in) :: block real (kind=RKIND), intent(in) :: dt @@ -937,7 +1063,9 @@ subroutine ocn_time_integrator_rk4_compute_thick_tends(block, dt, rkSubstepWeigh real (kind=RKIND), dimension(:), pointer :: sshCur real (kind=RKIND), dimension(:, :), pointer :: layerThicknessCur, normalVelocityCur + real (kind=RKIND), dimension(:, :), pointer :: vertVelocityNonhydro real (kind=RKIND), dimension(:, :), pointer :: normalVelocityProvis, highFreqThicknessProvis + real (kind=RKIND), dimension(:, :), pointer :: vertVelocityNonhydroProvis integer :: iEdge integer, pointer :: nEdges, nVertLevels @@ -962,6 +1090,7 @@ subroutine ocn_time_integrator_rk4_compute_thick_tends(block, dt, rkSubstepWeigh call mpas_pool_get_array(statePool, 'layerThickness', layerThicknessCur, 1) call mpas_pool_get_array(statePool, 'ssh', sshCur, 1) call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocityCur, 1) + call mpas_pool_get_array(provisStatePool, 'vertVelocityNonhydro', vertVelocityNonhydro, 1) call mpas_pool_get_array(provisStatePool, 'normalVelocity', normalVelocityProvis, 1) call mpas_pool_get_array(provisStatePool, 'highFreqThickness', highFreqThicknessProvis, 1) @@ -984,13 +1113,12 @@ subroutine ocn_time_integrator_rk4_compute_thick_tends(block, dt, rkSubstepWeigh if (associated(highFreqThicknessProvis)) then call ocn_vert_transport_velocity_top(verticalMeshPool, & layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & - sshCur, rkSubstepWeight, & - vertAleTransportTop, err, highFreqThicknessProvis) + sshCur, rkSubstepWeight, vertAleTransportTop, err, & + highFreqThicknessProvis) else call ocn_vert_transport_velocity_top(verticalMeshPool, & layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & - sshCur, rkSubstepWeight, & - vertAleTransportTop, err) + sshCur, rkSubstepWeight, vertAleTransportTop, err) endif call ocn_tend_thick(tendPool, forcingPool) @@ -1011,6 +1139,7 @@ subroutine ocn_time_integrator_rk4_compute_tracer_tends(block, dt, rkSubstepWeig real (kind=RKIND), dimension(:), pointer :: sshCur real (kind=RKIND), dimension(:, :), pointer :: layerThicknessCur, normalVelocityCur + real (kind=RKIND), dimension(:, :), pointer :: vertVelocityNonhydro real (kind=RKIND), dimension(:, :), pointer :: normalVelocityProvis, highFreqThicknessProvis logical, pointer :: config_filter_btr_mode @@ -1039,6 +1168,7 @@ subroutine ocn_time_integrator_rk4_compute_tracer_tends(block, dt, rkSubstepWeig call mpas_pool_get_array(statePool, 'layerThickness', layerThicknessCur, 1) call mpas_pool_get_array(statePool, 'ssh', sshCur, 1) call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocityCur, 1) + call mpas_pool_get_array(provisStatePool, 'vertVelocityNonhydro', vertVelocityNonhydro, 1) call mpas_pool_get_array(provisStatePool, 'normalVelocity', normalVelocityProvis, 1) call mpas_pool_get_array(provisStatePool, 'highFreqThickness', highFreqThicknessProvis, 1) @@ -1060,13 +1190,12 @@ subroutine ocn_time_integrator_rk4_compute_tracer_tends(block, dt, rkSubstepWeig if (associated(highFreqThicknessProvis)) then call ocn_vert_transport_velocity_top(verticalMeshPool, & layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & - sshCur, rkSubstepWeight, & - vertAleTransportTop, err, highFreqThicknessProvis) + sshCur, rkSubstepWeight, vertAleTransportTop, err, & + highFreqThicknessProvis) else call ocn_vert_transport_velocity_top(verticalMeshPool, & layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & - sshCur, rkSubstepWeight, & - vertAleTransportTop, err) + sshCur, rkSubstepWeight, vertAleTransportTop, err) endif if (config_filter_btr_mode) then @@ -1093,6 +1222,7 @@ subroutine ocn_time_integrator_rk4_compute_tends(domain, block, dt, rkWeight, dm real (kind=RKIND), dimension(:), pointer :: sshCur real (kind=RKIND), dimension(:, :), pointer :: layerThicknessCur, normalVelocityCur + real (kind=RKIND), dimension(:, :), pointer :: vertVelocityNonhydro real (kind=RKIND), dimension(:, :), pointer :: normalVelocityProvis, highFreqThicknessProvis logical, pointer :: config_filter_btr_mode @@ -1115,6 +1245,7 @@ subroutine ocn_time_integrator_rk4_compute_tends(domain, block, dt, rkWeight, dm call mpas_pool_get_array(statePool, 'layerThickness', layerThicknessCur, 1) call mpas_pool_get_array(statePool, 'ssh', sshCur, 1) call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocityCur, 1) + call mpas_pool_get_array(statePool, 'vertVelocityNonhydro', vertVelocityNonhydro, 1) call mpas_pool_get_array(provisStatePool, 'normalVelocity', normalVelocityProvis, 1) call mpas_pool_get_array(provisStatePool, 'highFreqThickness', highFreqThicknessProvis, 1) @@ -1123,13 +1254,12 @@ subroutine ocn_time_integrator_rk4_compute_tends(domain, block, dt, rkWeight, dm if (associated(highFreqThicknessProvis)) then call ocn_vert_transport_velocity_top(verticalMeshPool, & layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, & - sshCur, rkWeight, & - vertAleTransportTop, err, highFreqThicknessProvis) + sshCur, rkWeight, vertAleTransportTop, err, & + highFreqThicknessProvis) else call ocn_vert_transport_velocity_top(verticalMeshPool, & layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, & - sshCur, rkWeight, & - vertAleTransportTop, err) + sshCur, rkWeight, vertAleTransportTop, err) endif call ocn_tend_vel(domain, tendPool, provisStatePool, forcingPool, 1, dminfo, dt) @@ -1137,13 +1267,11 @@ subroutine ocn_time_integrator_rk4_compute_tends(domain, block, dt, rkWeight, dm if (associated(highFreqThicknessProvis)) then call ocn_vert_transport_velocity_top(verticalMeshPool, & layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & - sshCur, rkWeight, & - vertAleTransportTop, err, highFreqThicknessProvis) + sshCur, rkWeight, vertAleTransportTop, err, highFreqThicknessProvis) else call ocn_vert_transport_velocity_top(verticalMeshPool, & layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & - sshCur, rkWeight, & - vertAleTransportTop, err) + sshCur, rkWeight, vertAleTransportTop, err) endif call ocn_tend_thick(tendPool, forcingPool) @@ -1175,6 +1303,8 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkSubstepWeight, real (kind=RKIND), dimension(:, :), pointer :: normalVelocityCur, normalVelocityProvis, normalVelocityTend real (kind=RKIND), dimension(:, :), pointer :: layerThicknessCur, layerThicknessProvis, layerThicknessTend real (kind=RKIND), dimension(:, :), pointer :: lowFreqDivergenceCur, lowFreqDivergenceProvis, lowFreqDivergenceTend + real (kind=RKIND), dimension(:, :), pointer :: vertVelocityNonhydroCur, vertVelocityNonhydroTend + real (kind=RKIND), dimension(:, :), pointer :: vertVelocityNonhydroProvis real (kind=RKIND), dimension(:, :, :), pointer :: tracersGroupCur, tracersGroupProvis, tracersGroupTend @@ -1237,6 +1367,12 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkSubstepWeight, call mpas_pool_get_array(tracersPool, 'activeTracers', activeTracersCur, 1) + if(config_enable_nonhydrostatic_mode) then + call mpas_pool_get_array(statePool, 'vertVelocityNonhydro', vertVelocityNonhydroCur, 1) + call mpas_pool_get_array(tendPool, 'vertVelocityNonhydroTend', vertVelocityNonhydroTend) + call mpas_pool_get_array(provisStatePool, 'vertVelocityNonhydro', vertVelocityNonhydroProvis, 1) + end if + ! the loop below is on all vertlevels but we used edgemask to set ! velocities to zero below topo (esp. edges between topo and ocean) !$omp parallel @@ -1260,6 +1396,19 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkSubstepWeight, !$omp end do !$omp end parallel + if (config_enable_nonhydrostatic_mode) then + !$omp parallel + !$omp do schedule(runtime) private(k) + do iCell = 1, nCells + do k = 1, maxLevelCell(iCell) + vertVelocityNonhydroProvis(k, iCell) = vertVelocityNonhydroCur(k, iCell) + rkSubstepWeight & + * vertVelocityNonhydroTend(k, iCell) + end do + end do + !$omp end do + !$omp end parallel + end if + call mpas_pool_begin_iteration(tracersPool) do while ( mpas_pool_get_next_member(tracersPool, groupItr) ) if ( groupItr % memberType == MPAS_POOL_FIELD ) then @@ -1309,6 +1458,17 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkSubstepWeight, end do !$omp end do !$omp end parallel + + if (config_enable_nonhydrostatic_mode) then + !$omp parallel + !$omp do schedule(runtime) + do iCell = 1, nCells + vertVelocityNonhydroProvis(:,iCell) = vertVelocityNonhydroCur(:, iCell) + end do + !$omp end do + !$omp end parallel + end if + end if if (config_prescribe_thickness) then @@ -1433,6 +1593,8 @@ subroutine ocn_time_integrator_rk4_accumulate_update(block, rkWeight, err)!{{{ real (kind=RKIND), dimension(:, :), pointer :: layerThicknessNew, layerThicknessTend real (kind=RKIND), dimension(:, :), pointer :: highFreqThicknessNew, highFreqThicknessTend real (kind=RKIND), dimension(:, :), pointer :: lowFreqDivergenceNew, lowFreqDivergenceTend + real (kind=RKIND), dimension(:, :), pointer :: vertVelocityNonhydroNew + real (kind=RKIND), dimension(:, :), pointer :: vertVelocityNonhydroTend real (kind=RKIND), dimension(:, :, :), pointer :: tracersGroupNew, tracersGroupTend @@ -1470,6 +1632,11 @@ subroutine ocn_time_integrator_rk4_accumulate_update(block, rkWeight, err)!{{{ call mpas_pool_get_array(tendPool, 'normalVelocity', normalVelocityTend) call mpas_pool_get_array(tendPool, 'layerThickness', layerThicknessTend) + if (config_enable_nonhydrostatic_mode) then + call mpas_pool_get_array(tendPool, 'vertVelocityNonhydroTend', vertVelocityNonhydroTend) + call mpas_pool_get_array(statePool, 'vertVelocityNonhydro', vertVelocityNonhydroNew, 2) + end if + call mpas_pool_get_array(tendPool, 'highFreqThickness', highFreqThicknessTend) call mpas_pool_get_array(tendPool, 'lowFreqDivergence', lowFreqDivergenceTend) @@ -1499,6 +1666,16 @@ subroutine ocn_time_integrator_rk4_accumulate_update(block, rkWeight, err)!{{{ !$omp end do !$omp end parallel + if (config_enable_nonhydrostatic_mode) then + !$omp parallel + !$omp do schedule(runtime) + do iCell = 1, nCells + vertVelocityNonhydroNew(:, iCell) = vertVelocityNonhydroNew(:, iCell) + rkWeight * vertVelocityNonhydroTend(:, iCell) + end do + !$omp end do + !$omp end parallel + end if + call mpas_pool_begin_iteration(tracersPool) do while ( mpas_pool_get_next_member(tracersPool, groupItr) ) if ( groupItr % memberType == MPAS_POOL_FIELD ) then @@ -1548,8 +1725,9 @@ subroutine ocn_time_integrator_rk4_accumulate_update(block, rkWeight, err)!{{{ end subroutine ocn_time_integrator_rk4_accumulate_update!}}} - subroutine ocn_time_integrator_rk4_cleanup(block, dt, err)!{{{ + subroutine ocn_time_integrator_rk4_cleanup(block, dt, dminfo, err)!{{{ type (block_type), intent(in) :: block + type (dm_info), intent(in) :: dminfo real (kind=RKIND), intent(in) :: dt integer, intent(out) :: err @@ -1561,6 +1739,8 @@ subroutine ocn_time_integrator_rk4_cleanup(block, dt, err)!{{{ type (mpas_pool_type), pointer :: tracersPool real (kind=RKIND), dimension(:, :), pointer :: layerThicknessNew, normalVelocityNew + real (kind=RKIND), dimension(:, :), pointer :: vertVelocityNonhydroNew + real (kind=RKIND), dimension(:, :), pointer :: nonhydrostaticPressure real (kind=RKIND), dimension(:, :, :), pointer :: tracersGroupNew integer, dimension(:), pointer :: minLevelCell, maxLevelCell @@ -1613,6 +1793,10 @@ subroutine ocn_time_integrator_rk4_cleanup(block, dt, err)!{{{ call mpas_pool_get_array(forcingPool, 'landIceDraft', landIceDraft) endif + if (config_enable_nonhydrostatic_mode) then + call mpas_pool_get_array(statePool, 'vertVelocityNonhydro', vertVelocityNonhydroNew, 2) + end if + call mpas_pool_begin_iteration(tracersPool) do while ( mpas_pool_get_next_member(tracersPool, groupItr) ) if ( groupItr % memberType == MPAS_POOL_FIELD ) then @@ -1701,6 +1885,11 @@ subroutine ocn_time_integrator_rk4_cleanup(block, dt, err)!{{{ call ocn_vmix_implicit(dt, meshPool, statePool, forcingPool, scratchPool, err, 2) + if (config_enable_nonhydrostatic_mode) then + call ocn_nonhydrostatic_pressure_tend(normalVelocityNew, vertVelocityNonhydroNew, layerThicknessNew, & + layerThickEdgeMean, nhPressureCorrection, dminfo % comm, dt) + end if + end subroutine ocn_time_integrator_rk4_cleanup!}}} !*********************************************************************** diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F index a146560f23df..e39fabf587f6 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F @@ -45,6 +45,7 @@ module ocn_time_integration_si use ocn_gm use ocn_submesoscale_eddies + use ocn_compute_nonhydrostatic_pressure use ocn_equation_of_state use ocn_vmix use ocn_vertical_advection @@ -328,7 +329,11 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ highFreqThicknessNew, &! high frq thick at new time lowFreqDivergenceCur, &! low frq div at current time lowFreqDivergenceNew, &! low frq div at new time - tracerGroupIdealAgeMask ! mask for resetting surface ideal age + tracerGroupIdealAgeMask, &! mask for resetting surface ideal age + vertVelocityNonhydroNew, &! Vertical Velocity at new time + vertVelocityNonhydroCur, &! Vertical Velocity at current time + nonhydrostaticPressureNew, &! nhpressure at new time + nonhydrostaticPressureCur ! nh pressure at current time type (field1DReal), pointer :: & effectiveDensityField ! field pointer for halo update @@ -363,7 +368,8 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ normalVelocityTend, &! normal velocity tendency highFreqThicknessTend, &! thickness tendency for high-freq iter lowFreqDivergenceTend, &! thickness tendency for low freq - layerThicknessTend ! main thickness tendency + layerThicknessTend, &! main thickness tendency + vertVelocityNonhydroTend ! tendency of vertical velocity real (kind=RKIND), dimension(:,:,:), pointer :: & tracersGroupTend, &! tendencies for tracer groups @@ -504,7 +510,12 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ highFreqThicknessTend) call mpas_pool_get_array(tendPool, 'lowFreqDivergence', & lowFreqDivergenceTend) - + call mpas_pool_get_array(tendPool, 'vertVelocityNonhydroTend', & + vertVelocityNonhydroTend) + call mpas_pool_get_array(statePool, 'vertVelocityNonhydro', & + vertVelocityNonhydroNew, 2) + call mpas_pool_get_array(statePool, 'vertVelocityNonhydro', & + vertVelocityNonhydroCur, 1) call mpas_pool_get_array(tracersPool, 'activeTracers', activeTracersNew, 2) allocate(bottomDepthEdge(nEdgesAll+1)) @@ -642,6 +653,18 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ end if end do + if(config_enable_nonhydrostatic_mode) then + !$omp parallel + !$omp do schedule(runtime) private(k) + do iCell=1,nCellsAll + do k = 2, maxLevelCell(iCell) + vertVelocityNonhydroNew(k,iCell) = vertVelocityNonhydroCur(k,iCell) + end do + end do + !$omp end do + !$omp end parallel + end if + if (associated(highFreqThicknessNew)) then #ifdef MPAS_OPENACC !$acc parallel loop gang & @@ -825,13 +848,18 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ if (associated(highFreqThicknessNew)) then call ocn_vert_transport_velocity_top( & verticalMeshPool, layerThicknessCur, & - layerThickEdgeFlux, normalVelocityCur, sshCur, dt, & - vertAleTransportTop, err, highFreqThicknessNew) + layerThickEdgeFlux, normalVelocityCur, & + sshCur, dt, vertAleTransportTop, & + err, highFreqThicknessNew) + +#ifdef MPAS_OPENACC + !$acc exit data delete(highFreqThicknessNew) +#endif else call ocn_vert_transport_velocity_top( & verticalMeshPool, layerThicknessCur, & - layerThickEdgeFlux, normalVelocityCur, sshCur, dt, & - vertAleTransportTop, err) + layerThickEdgeFlux, normalVelocityCur, & + sshCur, dt, vertAleTransportTop, err) endif #ifdef MPAS_OPENACC @@ -2160,6 +2188,11 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ meshPool, swForcingPool, & dt, activeTracersOnly, 2) + if(config_enable_nonhydrostatic_mode) then + call ocn_tend_vvel(tendPool, statePool, & + forcingPool, meshPool, 2, dt) + end if + call mpas_pool_get_array(tracersPool, 'activeTracers', & tracersGroupCur, 1) call mpas_pool_get_array(tracersPool, 'activeTracers', & @@ -2240,6 +2273,12 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ call mpas_timer_stop("si halo tracers") call mpas_timer_start('si loop fini') + + if(config_enable_nonhydrostatic_mode) then + call mpas_pool_get_array(statePool, 'vertVelocityNonhydro', vertVelocityNonhydroCur, 1) + call mpas_pool_get_array(statePool, 'vertVelocityNonhydro', vertVelocityNonhydroNew, 2) + call mpas_pool_get_array(tendPool, 'vertVelocityNonhydroTend', vertVelocityNonhydroTend) + end if !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! If iterating, reset variables for next iteration @@ -2256,6 +2295,19 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ ! Only need T & S for earlier iterations, ! then all the tracers needed the last time through. + if(config_enable_nonhydrostatic_mode) then + !$omp parallel + !$omp do schedule(runtime) private(k, temp_h) + do iCell=1,nCells + do k=2,maxLevelCell(iCell) + temp_h = vertVelocityNonhydroCur(k,iCell) + dt * vertVelocityNonhydroTend(k,iCell) + vertVelocityNonhydroNew(k,iCell) = 0.5_RKIND*(vertVelocityNonhydroCur(k,iCell) + temp_h) + end do + end do + !$omp end do + !$omp end parallel + end if + !! Keeping this calc on the CPU for now #ifdef MPAS_OPENACC !$acc update host(layerThicknessNew) @@ -2413,6 +2465,18 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ !$omp end parallel #endif + if(config_enable_nonhydrostatic_mode) then + !$omp parallel + !$omp do schedule(runtime) private(k) + do iCell=1,nCells + do k=1,maxLevelCell(iCell) + vertVelocityNonhydroNew(k,iCell) = vertVelocityNonhydroCur(k,iCell) + dt*vertVelocityNonhydroTend(k,iCell) + end do + end do + !$omp end do + !$omp end parallel + end if + if (config_compute_active_tracer_budgets) then #ifdef MPAS_OPENACC @@ -2903,6 +2967,27 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ timeLevel=2) call mpas_timer_stop('si vmix halos normalVelFld') + if(config_enable_nonhydrostatic_mode) then + call mpas_pool_get_array(statePool, 'vertVelocityNonhydro', vertVelocityNonhydroNew, 2) + call mpas_pool_get_array(statePool, 'vertVelocityNonhydro', vertVelocityNonhydroCur, 1) + call mpas_pool_get_array(statePool, 'nonhydrostaticPressure', nonhydrostaticPressureNew, 2) + call mpas_pool_get_array(statePool, 'nonhydrostaticPressure', nonhydrostaticPressureCur, 1) + + dminfo = domain % dminfo + call ocn_nonhydrostatic_pressure_tend(normalVelocityNew, vertVelocityNonhydroNew, & + layerThicknessNew, layerThickEdgeMean, nhPressureCorrection, dminfo % comm, dt) + + call mpas_timer_start('se NHpressure tend halo pressure') + call mpas_dmpar_field_halo_exch(domain, 'nonhydrostaticPressure', timeLevel=2) + call mpas_dmpar_field_halo_exch(domain, 'nonhydrostaticPressure', timelevel=1) + call mpas_dmpar_field_halo_exch(domain, 'nhPressureCorrection') + call mpas_timer_stop('se NHpressure tend halo pressure') + + call ocn_nonhydrostatic_pressure_update_velocity(normalVelocityNew, vertVelocityNonhydroNew, & + layerThicknessNew, layerThickEdgeMean, nonhydrostaticPressureNew, nonhydrostaticPressureCur, & + nhPressureCorrection, dt) + end if + call mpas_pool_begin_iteration(tracersPool) do while ( mpas_pool_get_next_member(tracersPool, groupItr) ) if (groupItr%memberType == MPAS_POOL_FIELD) then diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F index ae171995ad80..9a251871e0b0 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F @@ -39,6 +39,7 @@ module ocn_time_integration_split use ocn_gm use ocn_submesoscale_eddies + use ocn_compute_nonhydrostatic_pressure use ocn_equation_of_state use ocn_vmix use ocn_vertical_advection @@ -136,6 +137,8 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ ! Local variables !----------------------------------------------------------------- + type (dm_info) :: dminfo + type (block_type), pointer :: & block ! structure with subdomain data @@ -218,7 +221,11 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ highFreqThicknessNew, &! high frq thick at new time lowFreqDivergenceCur, &! low frq div at current time lowFreqDivergenceNew, &! low frq div at new time - tracerGroupIdealAgeMask ! mask for resetting surface ideal age + tracerGroupIdealAgeMask, & ! mask for resetting surface ideal age + vertVelocityNonhydroCur, &! vertical velocity at current time + vertVelocityNonhydroNew, &! vertical velocity at new time + nonhydrostaticPressureNew, &! nonhydro pressure at new time + nonhydrostaticPressureCur ! nonhdyro pressure at current time type (field1DReal), pointer :: & effectiveDensityField ! field pointer for halo update @@ -250,7 +257,8 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ normalVelocityTend, &! normal velocity tendency highFreqThicknessTend, &! thickness tendency for high-freq iter lowFreqDivergenceTend, &! thickness tendency for low freq - layerThicknessTend ! main thickness tendency + layerThicknessTend, &! main thickness tendency + vertVelocityNonhydroTend ! vertical velocity tendency real (kind=RKIND), dimension(:,:,:), pointer :: & tracersGroupTend, &! tendencies for tracer groups @@ -328,6 +336,12 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ tracersTendPool) !*** Retrieve state variables at two time levels + if(config_enable_nonhydrostatic_mode) then + call mpas_pool_get_array(statePool, 'vertVelocityNonhydro', vertVelocityNonhydroCur, 1) + call mpas_pool_get_array(statePool, 'vertVelocityNonhydro', vertVelocityNonhydroNew, 2) + call mpas_pool_get_array(statePool, 'nonhydrostaticPressure', nonhydrostaticPressureNew, 2) + call mpas_pool_get_array(statePool, 'nonhydrostaticPressure', nonhydrostaticPressureCur, 1) + end if call mpas_pool_get_array(statePool, 'normalBaroclinicVelocity', & normalBaroclinicVelocityCur, 1) @@ -559,6 +573,18 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ #endif endif + if(config_enable_nonhydrostatic_mode) then + !$omp parallel + !$omp do schedule(runtime) private(k) + do iCell = 1, nCellsAll + do k = 2, maxLevelCell(iCell) + vertVelocityNonhydroNew(k,iCell) = vertVelocityNonhydroCur(k,iCell) + end do + end do + !$omp end do + !$omp end parallel + end if + call mpas_timer_stop("se prep") call mpas_pool_get_array(forcingPool, 'seaIcePressure', seaIcePressure) @@ -682,6 +708,11 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ call mpas_timer_start("se bcl vel") call mpas_timer_start('se bcl vel tend') + call mpas_pool_get_array(statePool, 'normalVelocity', & + normalVelocityCur, stage1_tend_time) + call mpas_pool_get_array(statePool, 'vertVelocityNonhydro', & + vertVelocityNonhydroCur, stage1_tend_time) + ! compute vertAleTransportTop. Use u (rather than & ! normalTransportVelocity) for momentum advection. ! Use the most recent time level available. @@ -689,13 +720,17 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ if (associated(highFreqThicknessNew)) then call ocn_vert_transport_velocity_top( & verticalMeshPool, layerThicknessCur, & - layerThickEdgeFlux, normalVelocityCur, sshCur, dt, & - vertAleTransportTop, err, highFreqThicknessNew) + layerThickEdgeFlux, normalVelocityCur, & + sshCur, dt, vertAleTransportTop, & + err, highFreqThicknessNew) +#ifdef MPAS_OPENACC + !$acc exit data delete(highFreqThicknessNew) +#endif else call ocn_vert_transport_velocity_top( & verticalMeshPool, layerThicknessCur, & - layerThickEdgeFlux, normalVelocityCur, sshCur, dt, & - vertAleTransportTop, err) + layerThickEdgeFlux, normalVelocityCur, & + sshCur, dt, vertAleTransportTop, err) endif #ifdef MPAS_OPENACC @@ -1668,16 +1703,17 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ #endif call ocn_vert_transport_velocity_top( & verticalMeshPool, layerThicknessCur, & - layerThickEdgeFlux, normalTransportVelocity, sshCur, & - dt, vertAleTransportTop, err, highFreqThicknessNew) + layerThickEdgeFlux, normalTransportVelocity, & + sshCur, dt, vertAleTransportTop, & + err, highFreqThicknessNew) #ifdef MPAS_OPENACC !$acc exit data delete(highFreqThicknessNew) #endif else call ocn_vert_transport_velocity_top( & verticalMeshPool, layerThicknessCur, & - layerThickEdgeFlux, normalTransportVelocity, sshCur, & - dt, vertAleTransportTop, err) + layerThickEdgeFlux, normalTransportVelocity, & + sshCur, dt, vertAleTransportTop, err) endif #ifdef MPAS_OPENACC !$acc exit data delete(layerThicknessCur, sshCur) @@ -1778,8 +1814,22 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ end if end do call mpas_timer_stop("se halo tracers") + if(config_enable_nonhydrostatic_mode) call ocn_tend_vvel(tendPool, statePool, & + forcingPool, meshPool, 2, dt) call mpas_timer_start('se loop fini') + call mpas_pool_get_array(tracersPool, 'activeTracers', & + tracersGroupCur, 1) + call mpas_pool_get_array(tracersPool, 'activeTracers', & + tracersGroupNew, 2) + if(config_enable_nonhydrostatic_mode) then + call mpas_pool_get_array(statePool, 'vertVelocityNonhydro', & + vertVelocityNonhydroCur, 1) + call mpas_pool_get_array(statePool, 'vertVelocityNonhydro', & + vertVelocityNonhydroNew, 2) + call mpas_pool_get_array(tendPool, 'vertVelocityNonhydroTend', & + vertVelocityNonhydroTend) + end if !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! If iterating, reset variables for next iteration @@ -1838,6 +1888,19 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ !$acc update device(activeTracersNew) #endif + if(config_enable_nonhydrostatic_mode) then + !$omp parallel + !$omp do schedule(runtime) private(k, temp_h) + do iCell = 1, nCells + do k = 2, maxLevelCell(iCell) + temp_h = vertVelocityNonhydroCur(k,iCell) + dt * vertVelocityNonhydroTend(k,iCell) + vertVelocityNonhydroNew(k,iCell) = 0.5_RKIND*(vertVelocityNonhydroCur(k,iCell) + temp_h) + end do + end do + !$omp end do + !$omp end parallel + end if + if (config_use_freq_filtered_thickness) then #ifdef MPAS_OPENACC @@ -1954,6 +2017,18 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ !$omp end parallel #endif + if(config_enable_nonhydrostatic_mode) then + !$omp parallel + !$omp do schedule(runtime) private(k) + do iCell = 1, nCells + do k = 1, maxLevelCell(iCell) + vertVelocityNonhydroNew(k,iCell) = vertVelocityNonhydroCur(k,iCell) + dt * vertVelocityNonhydroTend(k,iCell) + end do + end do + !$omp end do + !$omp end parallel + end if + if (config_compute_active_tracer_budgets) then #ifdef MPAS_OPENACC @@ -2430,6 +2505,25 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ ! and tracer fields. So this update is required to communicate ! the change due to implicit vertical mixing across the boundary. + if(config_enable_nonhydrostatic_mode) then + dminfo = domain % dminfo + call ocn_nonhydrostatic_pressure_tend(normalVelocityNew, & + vertVelocityNonhydroNew, layerThicknessNew, layerThickEdgeFlux, & + nhPressureCorrection, dminfo % comm, dt) + + call mpas_timer_start('se NH pressure tend halo pressure') + call mpas_dmpar_field_halo_exch(domain, 'nonhydrostaticPressure', timeLevel=2) + call mpas_dmpar_field_halo_exch(domain, 'nonhydrostaticPressure', timeLevel=1) + call mpas_dmpar_field_halo_exch(domain, 'nhPressureCorrection') + call mpas_timer_stop('se NH pressure tend halo pressure') + + call ocn_nonhydrostatic_pressure_update_velocity(normalVelocityNew, & + vertVelocityNonhydroNew, layerThicknessNew, layerThickEdgeFlux, & + nonhydrostaticPressureNew, nonhydrostaticPressureCur, & + nhPressureCorrection, dt) + end if + + call mpas_timer_start('se vmix halos') call mpas_timer_start('se vmix halos normalVelFld') @@ -2437,10 +2531,14 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ timeLevel=2) call mpas_timer_stop('se vmix halos normalVelFld') + call mpas_dmpar_field_halo_exch(domain, 'vertVelocityNonhydro', & + timeLevel=2) + call mpas_pool_begin_iteration(tracersPool) do while ( mpas_pool_get_next_member(tracersPool, groupItr) ) if (groupItr%memberType == MPAS_POOL_FIELD) then - + call mpas_dmpar_field_halo_exch(domain, & + groupItr%memberName, timeLevel=2) ! Reset iAge to zero where mask == 0 if (config_use_idealAgeTracers.and.trim(groupItr%memberName) == 'idealAgeTracers') then call mpas_pool_get_array(tracersPool, & diff --git a/components/mpas-ocean/src/ocean.cmake b/components/mpas-ocean/src/ocean.cmake index c1b71dcd5806..f8894aa6586d 100644 --- a/components/mpas-ocean/src/ocean.cmake +++ b/components/mpas-ocean/src/ocean.cmake @@ -8,10 +8,15 @@ list(APPEND INCLUDES "core_ocean/gotm/include") # check if lapack is linked find_package(LAPACK) find_package(BLAS) +find_package(PETSC) if(LAPACK_FOUND AND BLAS_FOUND) list(APPEND CPPDEFS "-DUSE_LAPACK") list(APPEND SLIBS "${LAPACK_LIBRARIES} ${BLAS_LIBRARIES}") endif() +if(PETSC_FOUND) + list(APPEND CPPDEFS "-DUSE_PETSC") + list(APPEND INCLUDES "${PETSC_PATH}/include") +endif() # driver (files live in E3SM) list(APPEND RAW_SOURCES @@ -104,6 +109,14 @@ list(APPEND RAW_SOURCES core_ocean/shared/mpas_ocn_time_varying_forcing.F core_ocean/shared/mpas_ocn_wetting_drying.F core_ocean/shared/mpas_ocn_vel_tidal_potential.F + core_ocean/shared/mpas_ocn_vvel_hmix_del2.F + core_ocean/shared/mpas_ocn_vvel_hmix_del4.F + core_ocean/shared/mpas_ocn_nonhydrostatic_pressure_solve.F + core_ocean/shared/mpas_ocn_vvel_pgrad.F + core_ocean/shared/mpas_ocn_vvel_coriolis.F + core_ocean/shared/mpas_ocn_vvel_advection.F + core_ocean/shared/mpas_ocn_vvel_horiz_advection.F + core_ocean/shared/mpas_ocn_vvel_vert_advection.F core_ocean/shared/mpas_ocn_stokes_drift.F ) diff --git a/components/mpas-ocean/src/shared/Makefile b/components/mpas-ocean/src/shared/Makefile index a9100a8bf9fe..ab69747fa3cf 100644 --- a/components/mpas-ocean/src/shared/Makefile +++ b/components/mpas-ocean/src/shared/Makefile @@ -27,6 +27,14 @@ OBJS = mpas_ocn_init_routines.o \ mpas_ocn_vel_pressure_grad.o \ mpas_ocn_vertical_remap.o \ mpas_ocn_vertical_regrid.o \ + mpas_ocn_vvel_vert_advection.o \ + mpas_ocn_vvel_horiz_advection.o \ + mpas_ocn_vvel_coriolis.o \ + mpas_ocn_vvel_hmix_del2.o \ + mpas_ocn_vvel_hmix_del4.o \ + mpas_ocn_vvel_advection.o \ + mpas_ocn_vvel_pgf.o \ + mpas_ocn_nonhydrostatic_pressure_solve.o \ mpas_ocn_vmix.o \ mpas_ocn_vmix_coefs_redi.o \ mpas_ocn_vmix_cvmix.o \ @@ -81,7 +89,7 @@ all: $(OBJS) mpas_ocn_init_routines.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o mpas_ocn_diagnostics.o mpas_ocn_diagnostics_variables.o mpas_ocn_gm.o mpas_ocn_submesoscale_eddies.o mpas_ocn_forcing.o mpas_ocn_surface_land_ice_fluxes.o -mpas_ocn_tendency.o: mpas_ocn_high_freq_thickness_hmix_del2.o mpas_ocn_tracer_surface_restoring.o mpas_ocn_thick_surface_flux.o mpas_ocn_tracer_short_wave_absorption.o mpas_ocn_tracer_advection.o mpas_ocn_tracer_hmix.o mpas_ocn_tracer_nonlocalflux.o mpas_ocn_surface_bulk_forcing.o mpas_ocn_surface_land_ice_fluxes.o mpas_ocn_tracer_surface_flux_to_tend.o mpas_ocn_tracer_interior_restoring.o mpas_ocn_tracer_exponential_decay.o mpas_ocn_tracer_ideal_age.o mpas_ocn_tracer_TTD.o mpas_ocn_vmix.o mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_frazil_forcing.o mpas_ocn_tidal_forcing.o mpas_ocn_tracer_ecosys.o mpas_ocn_tracer_DMS.o mpas_ocn_tracer_MacroMolecules.o mpas_ocn_tracer_CFC.o mpas_ocn_diagnostics.o mpas_ocn_wetting_drying.o mpas_ocn_vel_self_attraction_loading.o mpas_ocn_vel_tidal_potential.o mpas_ocn_mesh.o mpas_ocn_diagnostics_variables.o mpas_ocn_thick_hadv.o mpas_ocn_thick_vadv.o mpas_ocn_vel_hadv_coriolis.o mpas_ocn_vel_pressure_grad.o mpas_ocn_vel_vadv.o mpas_ocn_vel_hmix.o mpas_ocn_vel_forcing.o +mpas_ocn_tendency.o: mpas_ocn_high_freq_thickness_hmix_del2.o mpas_ocn_tracer_surface_restoring.o mpas_ocn_thick_surface_flux.o mpas_ocn_tracer_short_wave_absorption.o mpas_ocn_tracer_advection.o mpas_ocn_tracer_hmix.o mpas_ocn_tracer_nonlocalflux.o mpas_ocn_surface_bulk_forcing.o mpas_ocn_surface_land_ice_fluxes.o mpas_ocn_tracer_surface_flux_to_tend.o mpas_ocn_tracer_interior_restoring.o mpas_ocn_tracer_exponential_decay.o mpas_ocn_tracer_ideal_age.o mpas_ocn_tracer_TTD.o mpas_ocn_vmix.o mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_frazil_forcing.o mpas_ocn_tidal_forcing.o mpas_ocn_tracer_ecosys.o mpas_ocn_tracer_DMS.o mpas_ocn_tracer_MacroMolecules.o mpas_ocn_tracer_CFC.o mpas_ocn_diagnostics.o mpas_ocn_wetting_drying.o mpas_ocn_vel_self_attraction_loading.o mpas_ocn_vel_tidal_potential.o mpas_ocn_mesh.o mpas_ocn_diagnostics_variables.o mpas_ocn_thick_hadv.o mpas_ocn_thick_vadv.o mpas_ocn_vel_hadv_coriolis.o mpas_ocn_vel_pressure_grad.o mpas_ocn_vel_vadv.o mpas_ocn_vel_hmix.o mpas_ocn_vel_forcing.o mpas_ocn_vvel_pgrad.o mpas_ocn_vvel_advection.o mpas_ocn_vvel_coriolis.o mpas_ocn_vvel_hmix_del2.o mpas_ocn_diagnostics.o: mpas_ocn_thick_ale.o mpas_ocn_equation_of_state.o mpas_ocn_gm.o mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o mpas_ocn_diagnostics_variables.o mpas_ocn_surface_land_ice_fluxes.o mpas_ocn_vertical_advection.o mpas_ocn_submesoscale_eddies.o @@ -111,7 +119,7 @@ mpas_ocn_vel_vadv.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o mpas_ocn_vel_hmix.o: mpas_ocn_vel_hmix_del2.o mpas_ocn_vel_hmix_leith.o mpas_ocn_vel_hmix_del4.o mpas_ocn_constants.o mpas_ocn_config.o -mpas_ocn_vel_hmix_del2.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o +mpas_ocn_vel_hmix_del2.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o mpas_ocn_diagnostics_variables.o mpas_ocn_vel_hmix_leith.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o @@ -127,9 +135,25 @@ mpas_ocn_vel_forcing_explicit_bottom_drag.o: mpas_ocn_constants.o mpas_ocn_confi mpas_ocn_vel_hadv_coriolis.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o +mpas_ocn_vvel_coriolis.o: mpas_ocn_constants.o mpas_ocn_config.o + +mpas_ocn_vvel_pgf.o: mpas_ocn_diagnostics_variables.o mpas_ocn_constants.o mpas_ocn_config.o + +mpas_ocn_vvel_vert_advection.o: mpas_ocn_mesh.o mpas_ocn_tracer_advection_shared.o mpas_ocn_constants.o mpas_ocn_config.o + +mpas_ocn_vvel_horiz_advection.o: mpas_ocn_tracer_advection_shared.o mpas_ocn_mesh.o mpas_ocn_constants.o mpas_ocn_diagnostics_variables.o mpas_ocn_config.o + +mpas_ocn_vvel_hmix_del2.o: mpas_ocn_diagnostics_variables.o mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o + +mpas_ocn_vvel_hmix_del4.o: mpas_ocn_diagnostics_variables.o mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o + +mpas_ocn_vvel_advection.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_vvel_horiz_advection.o mpas_ocn_vvel_vert_advection.o + +mpas_ocn_nonhydrostatic_pressure_solve.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o mpas_ocn_diagnostics_variables.o + mpas_ocn_tracer_hmix.o: mpas_ocn_tracer_hmix_del2.o mpas_ocn_tracer_hmix_del4.o mpas_ocn_tracer_hmix_redi.o mpas_ocn_constants.o mpas_ocn_config.o -mpas_ocn_tracer_hmix_del2.o: mpas_ocn_constants.o mpas_ocn_config.o +mpas_ocn_tracer_hmix_del2.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o mpas_ocn_diagnostics_variables.o mpas_ocn_tracer_hmix_del4.o: mpas_ocn_constants.o mpas_ocn_config.o diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F index 4bdc554d09d6..6bb1bd3e833c 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F @@ -202,6 +202,8 @@ module ocn_diagnostics_variables real (kind=RKIND), dimension(:), pointer :: SIvec_w0,SIvec_w1 ,SIvec_wh0,SIvec_wh1,SIvec_y0 real (kind=RKIND), dimension(:), pointer :: SIvec_z0,SIvec_z1 ,SIvec_zh0,SIvec_zh1 + real (kind=RKIND), dimension(:,:), pointer :: nhPressureCorrection + !-------------------------------------------------------------------- ! ! Public member functions @@ -639,6 +641,8 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e call mpas_pool_get_array(diagnosticsPool, 'surfacePressure', & surfacePressure) + call mpas_pool_get_array(diagnosticsPool, 'nhPressureCorrection', nhPressureCorrection) + ! Semi-implicit Array Pointer retrievals if (trim(config_time_integrator) == 'split_implicit') then call mpas_pool_get_array(diagnosticsPool, 'SIvec_r0', SIvec_r0) @@ -779,7 +783,8 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e !$acc normalVelocitySurfaceLayer, & !$acc transportVelocityMeridional, & !$acc penetrativeTemperatureFluxOBL, & - !$acc normalizedRelativeVorticityCell & + !$acc normalizedRelativeVorticityCell, & + !$acc nhPressureCorrection & !$acc ) end if if ( trim(config_ocean_run_mode) == 'forward' ) then @@ -1021,7 +1026,8 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ !$acc normalVelocitySurfaceLayer, & !$acc transportVelocityMeridional, & !$acc penetrativeTemperatureFluxOBL, & - !$acc normalizedRelativeVorticityCell & + !$acc normalizedRelativeVorticityCell, & + !$acc nhPressureCorrection & !$acc ) end if if ( trim(config_ocean_run_mode) == 'forward' ) then @@ -1221,7 +1227,8 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ normalVelocitySurfaceLayer, & transportVelocityMeridional, & penetrativeTemperatureFluxOBL, & - normalizedRelativeVorticityCell) + normalizedRelativeVorticityCell, & + nhPressureCorrection) end if if ( trim(config_ocean_run_mode) == 'forward' ) then nullify(wettingVelocityFactor) diff --git a/components/mpas-ocean/src/shared/mpas_ocn_nonhydrostatic_pressure_solve.F b/components/mpas-ocean/src/shared/mpas_ocn_nonhydrostatic_pressure_solve.F new file mode 100644 index 000000000000..9eeabcf5efb6 --- /dev/null +++ b/components/mpas-ocean/src/shared/mpas_ocn_nonhydrostatic_pressure_solve.F @@ -0,0 +1,731 @@ +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! ocn_compute_nonhydrostatic_pressure +! +!> \brief MPAS ocean module to solve for the nonhydrostatic pressure correction +!> \author Luke Van Roekel +!> \date January 2021 +!> \details +! This module will derive the nonhydrostatic pressure correction +! uses petsc with the ILU preconditioner for solving +! +!----------------------------------------------------------------------- + +#ifdef USE_PETSC +#define PETSC_USE_FORTRAN_MODULES 1 +!#include "petsc/finclude/petscsys.h" +!#include "petsc/finclude/petscvec.h" +!#include "petsc/finclude/petscmat.h" +!#include "petsc/finclude/petscpc.h" +!#include "petsc/finclude/petscksp.h" +#include "petsc/finclude/petsc.h" +#endif + +module ocn_compute_nonhydrostatic_pressure + + use mpas_derived_types + use mpas_pool_routines + use mpas_timer + use mpas_constants + use mpas_log + + use ocn_constants + use ocn_config + use ocn_mesh + use ocn_diagnostics_variables + +#ifdef USE_PETSC + !PETSC includes for the matrix solve + use petscsys + use petscvec + use petscmat + use petscpc + use petscksp + use petsc +#endif + implicit none + private + save + + !-------------------------------------------------------------------- + ! + ! Public member functions + ! + !-------------------------------------------------------------------- + + public :: ocn_nonhydrostatic_pressure_tend, & + ocn_nonhydrostatic_solver_init, & + ocn_nonhydrostatic_solver_finalize, & + ocn_nonhydrostatic_solver, & + ocn_nonhydrostatic_pressure_update_velocity + + !-------------------------------------------------------------------- + ! + ! Private module variables + ! + !-------------------------------------------------------------------- + + double precision :: norm ! norm of solution error +#ifdef USE_PETSC + PetscInt :: i, j, II, JJ, m, n, its, Istart, Iend + PetscInt :: row, col, ione, globalM, first, last, firstVal + PetscInt :: overlap,nlocal + PetscErrorCode :: petsc_err, ierr2 + PetscMPIInt :: petsc_rank, petsc_size + PetscBool :: flag + PetscScalar :: v, one, neg_one + Vec :: x, b, u, negOneVec ! approximate solution, right hand side vector, exact solution vector + Mat :: Amat ! Matrix that defines the system + Mat :: Pmat ! Preconditioner Matrix + KSP :: ksp ! krylov subspace method context + KSP,allocatable,dimension(:) :: subksp + PetscRandom :: rctx ! random number generator + PC :: pc, subpc ! preconditioner context + PCType :: ptype ! flag to set type of preconditioner MAYBE make this a NL option? + PetscReal :: tol + + !Stuff for the new vecscatter + IS :: origIS, destIS !PETSc index sets + VecScatter :: scatter ! PETSc vec scatter context computed once and reused + Vec :: local_x ! local value of things to grab from global vector + PetscInt,allocatable,dimension(:) :: destIndex, origIndex + PetscScalar,pointer :: values(:) +#endif + + real(kind=RKIND) :: surfaceBoundaryConditionTop, surfaceBoundaryConditionBottom + + contains + + subroutine ocn_nonhydrostatic_pressure_tend(normalVelocityNew, & + vertVelocityNonhydroNew, layerThickness, layerThicknessEdge, & + nhPressureCorrection, mpi_comm, dt) + + real(kind=RKIND), dimension(:,:), intent(in) :: normalVelocityNew, vertVelocityNonhydroNew, & + layerThickness, layerThicknessEdge + + real(kind=RKIND), intent(in) :: dt + + real(kind=RKIND), dimension(:,:), intent(inout) :: nhPressureCorrection + + integer, intent(in) :: mpi_comm + + logical :: reuse_preconditioner + + reuse_preconditioner = .false. + + call ocn_nonhydrostatic_solver(normalVelocityNew, & + vertVelocityNonhydroNew, layerThickness, layerThicknessEdge, mpi_comm, & + dt, nhPressureCorrection) + + end subroutine ocn_nonhydrostatic_pressure_tend + + !-------------------------------------------------------------------- + ! + ! routine ocn_nonhydrostatic_pressure_update_velocity + ! + !> \brief Update the normal and vertical velocity + !> \author Luke Van Roekel + !> \date February 2021 + !> \details + !> creates KSP solver context, gets preconditioner and solves + ! + !-------------------------------------------------------------------- + + subroutine ocn_nonhydrostatic_pressure_update_velocity(normalVelocityNew, & + vertVelocityNonhydroNew, layerThickness, layerThicknessEdge, & + nonhydrostaticPressure, nonhydrostaticPressureOld, & + nhPressureCorrection, dt) + + real(KIND=RKIND),intent(in) :: dt + + real(KIND=RKIND),dimension(:,:), intent(inout) :: & + normalVelocityNew, & + vertVelocityNonhydroNew + + real(KIND=RKIND),dimension(:,:),intent(in) :: & + layerThicknessEdge, & + layerThickness + + real(KIND=RKIND),dimension(:,:),intent(inout) :: & + nonhydrostaticPressure, & + nonhydrostaticPressureOld, & + nhPressureCorrection + + integer :: cell1, cell2, nEdges, iCell, iEdge, nCells, i, k + + real(kind=RKIND), dimension(nVertLevels) :: div_hu + + real(kind=RKIND) :: invAreaCell, r_tmp + + nEdges = nEdgesHalo(2) + + call mpas_timer_start('NH solve - update velocity',.false.) + ! loop over edges to update normalVelocity + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(k, cell1, cell2) + do iEdge=1,nEdges + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + do k = 1,maxLevelEdgeTop(iEdge) + normalVelocityNew(k,iEdge) = normalVelocityNew(k,iEdge) - dt * & + (nhPressureCorrection(k,cell2) - nhPressureCorrection(k,cell1)) / & + dcEdge(iEdge)*edgeMask(k,iEdge) + end do + end do + !$omp end do + !$omp end parallel + + nCells = nCellsHalo( 1 ) + + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(i, r_tmp, iEdge, invAreaCell, k, div_hu) + do iCell = 1, nCells + div_hu(:) = 0.0_RKIND + invAreaCell = 1.0_RKIND / areaCell(iCell) + do i = 1,nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + do k = 1,maxLevelEdgeTop(iEdge) + r_tmp = dvEdge(iEdge)*normalVelocityNew(k,iEdge)*invAreaCell + div_hu(k) = div_hu(k) - layerThicknessEdge(k,iEdge)*edgeSignOnCell(i,iCell)*r_tmp + end do + end do + + vertVelocityNonhydroNew(maxLevelCell(iCell) + 1, iCell) = 0.0_RKIND + do k=maxLevelCell(iCell),minLevelCell(iCell),-1 + vertVelocityNonhydroNew(k,iCell) = vertVelocityNonhydroNew(k+1,iCell) - div_hu(k) + end do + end do + !$omp end do + !$omp end parallel + + !Update the nonhydrostatic pressure to be second order accurate + + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(k) + do iCell=1,nCells + do k = minLevelCell(iCell),maxLevelCell(iCell) + nonhydrostaticPressure(k,iCell) = nonhydrostaticPressureOld(k,iCell) + & + nhPressureCorrection(k,iCell) + end do + end do + !$omp end do + !$omp end parallel + + call mpas_timer_stop('NH solve - update velocity') + + end subroutine ocn_nonhydrostatic_pressure_update_velocity + + !-------------------------------------------------------------------- + ! + ! routine ocn_nonhydrostatic_solver + ! + !> \brief Create the matrix and right-hand-side and then solve the + !> linear system + !> \author Sara Calandrini, Darren Engwirda, Luke Van Roekel + !> \date November 2021 + !> \details + !> creates KSP solver context, gets preconditioner and solves + ! + !-------------------------------------------------------------------- + + subroutine ocn_nonhydrostatic_solver(normalVelocityNew, & + vertVelocityNonhydroNew, layerThickness, layerThicknessEdge, mpi_comm, & + dt, nhPressureCorrection) + + implicit none + + real(kind=RKIND), dimension(:,:), intent(in) :: normalVelocityNew, & + vertVelocityNonhydroNew + real(kind=RKIND), dimension(:,:), intent(in) :: layerThickness, & + layerThicknessEdge + real(kind=RKIND), intent(in) :: dt + real(kind=RKIND), dimension(:,:), intent(inout) :: nhPressureCorrection + integer, intent(in) :: mpi_comm + + !-------------------------------------------------------------------- + ! + ! subroutine variables + ! + !-------------------------------------------------------------------- + + integer :: k, iCell, jCell, iEdge, jEdge, iCell3d, jCell3d + integer :: nnzRowMax, vecSpot + real(kind=RKIND) :: xEdge, thicknessUpper, thicknessLower + real(kind=RKIND), dimension(nVertLevels) :: div + +#ifdef USE_PETSC + PetscScalar :: insertVal, rhsSum + PetscInt :: itNum + + call mpas_timer_start('NH solve - matrix',.false.) + + ! form matrix: del^2(q) = 1/dt * div(u*) + + !nnzRowMax = 9 * 2 + 4 ! do this better + nnzRowMax = 1 + 2 + maxval(nEdgesOnCell) + + call MatZeroEntries(Amat, petsc_err) + + do iCell = 1, nCellsOwned + + ! x-part of del^2 operator + + do jEdge = 1, nEdgesOnCell(iCell) + + iEdge = edgesOnCell(jEdge,iCell) + jCell = cellsOnCell(jEdge,iCell) + + xEdge = dvEdge(iEdge) / dcEdge(iEdge) + + do k = minLevelCell(iCell), maxLevelCell(iCell) + if (boundaryEdge(k,iEdge) .ge. 1) cycle + + iCell3d = (indexToCellId(iCell) - 1) * nVertLevels + k - 1 + jCell3d = (indexToCellId(jCell) - 1) * nVertLevels + k - 1 + + insertVal = +1 * layerThicknessEdge(k,iEdge) * xEdge + call MatSetValues(Amat, 1, iCell3d, 1, jCell3d, insertVal, ADD_VALUES, petsc_err) + insertVal = -1 * layerThicknessEdge(k,iEdge) * xEdge + call MatSetValues(Amat, 1, iCell3d, 1, iCell3d, insertVal, ADD_VALUES, petsc_err) + end do + + end do + + ! z-part of del^2 operator + + ! first layer upper q = 0 + ! first layer lower + k = minLevelCell(iCell) !1 + + thicknessUpper = 1.0E-15_RKIND + 0.5_RKIND * ( & + layerThickness(k,iCell) ) + thicknessLower = 1.0E-15_RKIND + 0.5_RKIND * ( & + layerThickness(k,iCell) + layerThickness(k+1,iCell)) + + iCell3d = (indexToCellId(iCell) - 1) * nVertLevels + k - 1 + + insertVal = surfaceBoundaryConditionTop*(-1 * areaCell(iCell) / thicknessUpper) + call MatSetValues(Amat, 1, iCell3d, 1, iCell3d, insertVal, ADD_VALUES, petsc_err) + + iCell3d = (indexToCellId(iCell) - 1) * nVertLevels + k - 1 + jCell3d = (indexToCellId(iCell) - 1) * nVertLevels + (k+1) - 1 + + insertVal = +1 * areaCell(iCell) / thicknessLower + call MatSetValues(Amat, 1, iCell3d, 1, jCell3d, insertVal, ADD_VALUES, petsc_err) + insertVal = -1 * areaCell(iCell) / thicknessLower + call MatSetValues(Amat, 1, iCell3d, 1, iCell3d, insertVal, ADD_VALUES, petsc_err) + + do k = 2, maxLevelCell(iCell) - 1 + + thicknessUpper = 1.0E-15_RKIND + 0.5_RKIND * ( & + layerThickness(k,iCell) + layerThickness(k-1,iCell)) + thicknessLower = 1.0E-15_RKIND + 0.5_RKIND * ( & + layerThickness(k,iCell) + layerThickness(k+1,iCell)) + + ! upper + iCell3d = (indexToCellId(iCell) - 1) * nVertLevels + k - 1 + jCell3d = (indexToCellId(iCell) - 1) * nVertLevels + (k-1) - 1 + + insertVal = +1 * areaCell(iCell) / thicknessUpper + call MatSetValues(Amat, 1, iCell3d, 1, jCell3d, insertVal, ADD_VALUES, petsc_err) + insertVal = -1 * areaCell(iCell) / thicknessUpper + call MatSetValues(Amat, 1, iCell3d, 1, iCell3d, insertVal, ADD_VALUES, petsc_err) + + ! lower + iCell3d = (indexToCellId(iCell) - 1) * nVertLevels + k - 1 + jCell3d = (indexToCellId(iCell) - 1) * nVertLevels + (k+1) - 1 + + insertVal = +1 * areaCell(iCell) / thicknessLower + call MatSetValues(Amat, 1, iCell3d, 1, jCell3d, insertVal, ADD_VALUES, petsc_err) + insertVal = -1 * areaCell(iCell) / thicknessLower + call MatSetValues(Amat, 1, iCell3d, 1, iCell3d, insertVal, ADD_VALUES, petsc_err) + + end do + + ! final layer upper + ! final layer lower dq/dz = 0 + k = maxLevelCell(iCell) + + thicknessUpper = 1.0E-15_RKIND + 0.5_RKIND * ( & + layerThickness(k,iCell) + layerThickness(k-1,iCell)) + thicknessLower = 1.0E-15_RKIND + 0.5_RKIND * ( & + layerThickness(k,iCell) ) + + iCell3d = (indexToCellId(iCell) - 1) * nVertLevels + k - 1 + + insertVal = surfaceBoundaryConditionBottom*(-1 * areaCell(iCell) / thicknessLower) + call MatSetValues(Amat, 1, iCell3d, 1, iCell3d, insertVal, ADD_VALUES, petsc_err) + + iCell3d = (indexToCellId(iCell) - 1) * nVertLevels + k - 1 + jCell3d = (indexToCellId(iCell) - 1) * nVertLevels + (k-1) - 1 + + insertVal = +1 * areaCell(iCell) / thicknessUpper + call MatSetValues(Amat, 1, iCell3d, 1, jCell3d, insertVal, ADD_VALUES, petsc_err) + insertVal = -1 * areaCell(iCell) / thicknessUpper + call MatSetValues(Amat, 1, iCell3d, 1, iCell3d, insertVal, ADD_VALUES, petsc_err) + + ! for non-flat bathymetry + do k = maxLevelCell(iCell)+1, nVertLevels + iCell3D = (indexToCellId(iCell) - 1) * nVertLevels + k - 1 + insertVal = 1.0_RKIND + call MatSetValues(Amat, 1, iCell3d, 1, iCell3d, insertVal, ADD_VALUES, petsc_err) + end do + + end do + + call MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY,petsc_err) + call MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY,petsc_err) + + call mpas_timer_stop('NH solve - matrix') + + ! form vector: del^2(q) = 1/dt * div(u*) + + call mpas_timer_start('NH solve - rhs',.false.) + + call vecZeroEntries(b,petsc_err) + + do iCell = 1, nCellsOwned + + div(:) = 0.0_RKIND + + ! x-part of 1/dt * div(u*) + + do jEdge = 1, nEdgesOnCell(iCell) + + iEdge = edgesOnCell(jEdge,iCell) + + do k = minLevelCell(iCell), maxLevelCell(iCell) + + div(k) = div(k) - & + edgeSignOnCell(jEdge,iCell) * dvEdge(iEdge) * & + normalVelocityNew(k,iEdge) * layerThicknessEdge(k,iEdge) + + end do + + end do + + ! z-part of 1/dt * div(u*) + + do k = minLevelCell(iCell), maxLevelCell(iCell) + + div(k) = div(k) + areaCell(iCell) * ( & + vertVelocityNonhydroNew(k,iCell) - vertVelocityNonhydroNew(k+1,iCell)) + + end do + + ! copy into full rhs vector + + do k = minLevelCell(iCell), maxLevelCell(iCell) + iCell3d = (indexToCellId(iCell) - 1) * nVertLevels + k - 1 + insertVal = div(k) / dt + call VecSetValues(b, 1, iCell3d, insertVal, INSERT_VALUES, petsc_err) + end do + + end do + + call VecAssemblyBegin(b,petsc_err) + call VecAssemblyEnd(b,petsc_err) + + if (config_nonhydrostatic_remove_rhs_mean) then + call VecSum(b, rhsSum, petsc_err) + rhsSum = rhsSum / globalM + call VecAXPY(b, rhsSum, negOneVec, petsc_err) + end if + + call mpas_timer_stop('NH solve - rhs') + + ! PCG solver: del^2(q) = 1/dt * div(u*) + + call mpas_timer_start('NH solve - solver',.false.) + + call KSPCreate(mpi_comm, ksp, petsc_err) + call KSPSetType(ksp, config_nonhydrostatic_solver_type, petsc_err) + call KSPSetOperators(ksp,Amat,Amat,petsc_err) + call KSPGetPC(ksp,pc,petsc_err) + call PCSetType(pc,config_nonhydrostatic_preconditioner,petsc_err) + call KSPSetFromOptions(ksp,petsc_err) + call KSPSetTolerances(ksp,config_petsc_rtol,config_petsc_atol,1.0E4_RKIND,config_petsc_maxit,petsc_err) + call KSPSolve(ksp, b, x, petsc_err) + call KSPGetIterationNumber(ksp, itNum, petsc_err) + print *, 'iteration count = ',itNum + call KSPSetInitialGuessNonzero(ksp,PETSC_TRUE,petsc_err) + call KSPDestroy(ksp,ierr2) + + !now update the nonhydrostaticPressure array unpack into k,iCell + call VecScatterBegin(scatter, x, local_x, INSERT_VALUES, SCATTER_FORWARD, petsc_err) + call VecScatterEnd(scatter, x, local_x, INSERT_VALUES, SCATTER_FORWARD, petsc_err) + call VecGetArrayReadF90(local_x, values, petsc_err) + + !$omp parallel + !$omp do schedule(runtime) private(k,vecSpot) + do iCell=1,nCellsOwned + do k=minLevelCell(iCell),maxLevelCell(iCell) + vecSpot = (iCell - 1) * nVertLevels + k - 1 + nhPressureCorrection(k,iCell) = values(vecSpot+1) + end do + end do + !print*, nhPressureCorrection(10,100) + !$omp end do + !$omp end parallel + + call mpas_timer_stop('NH solve - solver') +#endif + + end subroutine ocn_nonhydrostatic_solver + + + !-------------------------------------------------------------------- + ! €I€I€I€I€I€I + ! routine ocn_nonhydrostatic_solver_init + ! + !> \brief Initializes matrices€I€I€I€I€I€I and vectors for PETSC solve + !> \author Luke Van Roekel + !> \date February 2021 + !> \details + !> Initializese the matrix for the nonhydrostatic solve, automatically + !> determines what is owned by processor. Does not fill values. + ! + !-------------------------------------------------------------------- + + subroutine ocn_nonhydrostatic_solver_init(mpi_comm, ierr) + + integer, intent(in) :: mpi_comm + integer, intent(inout) :: ierr + integer :: commsize, globalCells + real(kind=RKIND) :: decimals, version +#ifdef USE_PETSC + PetscInt :: major, minor, subminor, release, locRow, locCol, locals(nCellsOwned) + PetscReal :: val + ISLocalToGlobalMapping :: mapping +#endif + integer :: k, iCell, iter, nnzRowMax + + nnzRowMax = 1 + 2 + maxval(nEdgesOnCell) + +#ifndef USE_PETSC + call mpas_log_write("nonhydrostatic solver requires PETSC library compiled and linked", & + mpas_log_crit) +#endif + + if(config_nonhydrostatic_solve_surface_boundary_condition == 'noGradient') then + surfaceBoundaryConditionTop = 0.0_RKIND + surfaceBoundaryConditionBottom = 0.0_RKIND + elseif(config_nonhydrostatic_solve_surface_boundary_condition == 'noPressure') then + surfaceBoundaryConditionTop = 1.0_RKIND + surfaceBoundaryConditionBottom = 1.0_RKIND + elseif(config_nonhydrostatic_solve_surface_boundary_condition == 'pressureTopGradientBottom') then + surfaceBoundaryConditionTop = 1.0_RKIND + surfaceBoundaryConditionBottom = 0.0_RKIND + else + surfaceBoundaryConditionTop = 0.0_RKIND + surfaceBoundaryConditionBottom = 1.0_RKIND + endif + + call MPI_comm_size(mpi_comm, commsize, ierr) + +#ifdef USE_PETSC + call PetscInitialize(PETSC_NULL_CHARACTER,petsc_err) + if(petsc_err .ne. 0) then + call mpas_log_write("error: petsc initialize failed, error code = $i", mpas_log_crit, & + intargs=(/petsc_err/)) + ierr = petsc_err + endif + + call PetscGetVersionNumber(major, minor, subminor, release, petsc_err) + if(petsc_err .ne. 0) then + call mpas_log_write("ERROR: Petsc GetVersionNumber failed, error code = $i", MPAS_LOG_CRIT, & + intArgs=(/petsc_err/)) + ierr = petsc_err + end if + + !Solver routines require version to be greater than 3.5 + decimals = floor(log10(float(minor))) + 1 + version = float(major) + float(minor) / 10.0_RKIND**decimals + if(version < 3.05) then + call mpas_log_write("ERROR: Petsc Version required to be greater than 3.5, you are using $r", & + MPAS_LOG_CRIT, realArgs=(/version/)) + ierr = petsc_err + end if + + m = nCellsOwned * nVertLevels + call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-m',m,flag,petsc_err) + if(petsc_err .ne. 0) then + call mpas_log_write("ERROR: Petsc OptionsGetInt failed, error code = $i", MPAS_LOG_CRIT, & + intArgs=(/petsc_err/)) + ierr = petsc_err + endif + + call MatCreate(mpi_comm, Amat, petsc_err) + if(petsc_err .ne. 0) then + call mpas_log_write("ERROR: Petsc MatCreate failed, error code = $i", MPAS_LOG_CRIT, & + intArgs=(/petsc_err/)) + ierr = petsc_err + endif + + call MatSetSizes(Amat, m, m, PETSC_DETERMINE, PETSC_DETERMINE, petsc_err) + if(petsc_err .ne. 0) then + call mpas_log_write("ERROR: Petsc MatSetSizes failed, error code = $i", MPAS_LOG_CRIT, & + intArgs=(/petsc_err/)) + ierr = petsc_err + endif + + call MatSetType( Amat, MATAIJ, petsc_err ) + if( commsize == 1) then + call MatSetType( Amat, MATAIJ, petsc_err ) + call MatSetFromOptions(Amat,petsc_err) + call MatSeqAIJSetPreallocation(Amat, nnzRowMax, PETSC_NULL_INTEGER, petsc_err) + else + call MatSetType( Amat, MATMPIAIJ, petsc_err ) + call MatSetFromOptions(Amat,petsc_err) + call MatMPIAIJSetPreallocation(Amat, nnzRowMax, PETSC_NULL_INTEGER, nnzRowMax, PETSC_NULL_INTEGER, petsc_err) + end if + + call MatGetOwnershipRange(Amat,Istart,Iend,petsc_err) + if(petsc_err .ne. 0) then + call mpas_log_write("ERROR: Petsc MatGetOwnershipRange failed, error_code = $i", MPAS_LOG_CRIT, & + intArgs=(/petsc_err/)) + ierr = petsc_err + endif + + call MatGetSize(Amat, globalM, globalM, petsc_err) + if(petsc_err .ne. 0) then + call mpas_log_write("ERROR: Petsc MatGetSize failed, error_code = $i", MPAS_LOG_CRIT, & + intArgs=(/petsc_err/)) + ierr = petsc_err + endif + + call MatGetLocalSize(Amat,locRow,locCol,petsc_err) + + print*, m, globalM, locRow, locCol + + call VecCreate(mpi_comm, b, petsc_err) + if(petsc_err .ne. 0) then + call mpas_log_write("ERROR: Petsc VecCreate failed, error_code = $i", MPAS_LOG_CRIT, & + intArgs=(/petsc_err/)) + ierr = petsc_err + endif + + call VecSetSizes(b, m, globalM, petsc_err) + if ( commsize == 1 ) then + call VecSetType(b, VECSEQ, petsc_err) + else + call VecSetType(b, VECMPI, petsc_err) + end if + + if(petsc_err .ne. 0) then + call mpas_log_write("ERROR: Petsc VecSetSizes failed, error_code = $i", MPAS_LOG_CRIT, & + intArgs=(/petsc_err/)) + ierr = petsc_err + endif + + call VecDuplicate(b, x, petsc_err) + if(petsc_err .ne. 0) then + call mpas_log_write("ERROR: Petsc VecCreate Duplicate, error_code = $i", MPAS_LOG_CRIT, & + intArgs=(/petsc_err/)) + ierr = petsc_err + endif + + call VecSetSizes(x, m, globalM, petsc_err) + if(petsc_err .ne. 0) then + call mpas_log_write("ERROR: Petsc VecSetSizes failed, error_code = $i", MPAS_LOG_CRIT, & + intArgs=(/petsc_err/)) + ierr = petsc_err + endif + + call VecGetOwnershipRange(x,first,last,petsc_err) + + call VecDuplicate(b, negOneVec, petsc_err) + if(petsc_err .ne. 0) then + call mpas_log_write("ERROR: Petsc VecCreate Duplicate, error_code = $i", MPAS_LOG_CRIT, & + intArgs=(/petsc_err/)) + ierr = petsc_err + endif + + call VecSetSizes(negOneVec, m, globalM, petsc_err) + if(petsc_err .ne. 0) then + call mpas_log_write("ERROR: Petsc VecSetSizes failed, error_code = $i", MPAS_LOG_CRIT, & + intArgs=(/petsc_err/)) + ierr = petsc_err + endif + + call VecGetOwnershipRange(negOneVec,first,last,petsc_err) + + call VecCreate(mpi_comm, local_x, petsc_err) + if(petsc_err .ne. 0) then + call mpas_log_write("ERROR: Petsc VecCreate failed, error_code = $i", MPAS_LOG_CRIT, & + intArgs=(/petsc_err/)) + ierr = petsc_err + endif + + call VecSetSizes(local_x, m, PETSC_DETERMINE, petsc_err) + if ( commsize == 1 ) then + call VecSetType(local_x, VECSEQ, petsc_err) + else + call VecSetType(local_x, VECMPI, petsc_err) + end if + + allocate(destIndex(m), origIndex(m)) + ! Try a better way with Index Sets and Vec Scatter for some communication + ! later move this first part to init as we don't have to do it every time FIXME! + iter = 1 + do iCell = 1,nCellsOwned + do k=1,nVertLevels + destIndex(iter) = (iCell-1)*nVertLevels + k - 1 + first + origIndex(iter) = (indexToCellID(iCell)-1)*nVertLevels + k - 1 + call VecSetValues(negOneVec, 1, origIndex(iter), -1.0_RKIND, INSERT_VALUES, petsc_err) + iter = iter + 1 + end do + end do + + call VecAssemblyBegin(negOneVec, petsc_err) + call VecAssemblyEnd(negOneVec, petsc_err) + + call ISCreateGeneral(mpi_comm, m, origIndex, PETSC_COPY_VALUES, origIS, petsc_err) + call ISCreateGeneral(mpi_comm, m, destIndex, PETSC_COPY_VALUES, destIS, petsc_err) + call VecScatterCreate(x, origIS, local_x, destIS, scatter, petsc_err) + + deallocate(destIndex, origIndex) + + call ISDestroy(origIS, petsc_err) + call ISDestroy(destIS, petsc_err) + +#endif + + end subroutine ocn_nonhydrostatic_solver_init + + !-------------------------------------------------------------------- + ! + ! routine ocn_nonhydrostatic_solver_finalize + ! + !> \brief Destroys matrices and vectors for PETSC solve + !> \author Luke Van Roekel + !> \date February 2021 + !> \details + !> Destroys matrices and frees memory, finalizes Petsc + ! + !-------------------------------------------------------------------- + + subroutine ocn_nonhydrostatic_solver_finalize() + +#ifdef USE_PETSC + !stuff for PETSC_finalize from the scatter routine + call VecScatterDestroy(scatter, petsc_err) + + call MatDestroy(Amat,petsc_err) + call VecDestroy(b,petsc_err) + call VecDestroy(x,petsc_err) + call VecDestroy(negOneVec,petsc_err) + call VecDestroy(local_x,petsc_err) + + !LPV FIXME : for some reason this causes error "Corrupt argument" + !disable for now + !call PetscFinalize(petsc_err) +#endif + + end subroutine ocn_nonhydrostatic_solver_finalize + +end module ocn_compute_nonhydrostatic_pressure + diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tendency.F b/components/mpas-ocean/src/shared/mpas_ocn_tendency.F index 2341ea5c28da..31adcba0d70f 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_tendency.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_tendency.F @@ -28,7 +28,6 @@ module ocn_tendency use mpas_threading use mpas_dmpar use ocn_diagnostics - use ocn_diagnostics_variables use ocn_constants use ocn_config use ocn_mesh @@ -37,6 +36,7 @@ module ocn_tendency use ocn_surface_land_ice_fluxes use ocn_frazil_forcing use ocn_tidal_forcing + use ocn_diagnostics_variables use ocn_tracer_hmix use ocn_high_freq_thickness_hmix_del2 @@ -68,6 +68,14 @@ module ocn_tendency use ocn_vel_tidal_potential use ocn_vel_self_attraction_loading + use ocn_vvel_pgrad + use ocn_vvel_coriolis + use ocn_vvel_advection + use ocn_vvel_hmix_del2 + use ocn_vvel_hmix_del4 + + use mpas_vector_reconstruction + implicit none private save @@ -88,7 +96,8 @@ module ocn_tendency ocn_tend_vel, & ocn_tend_tracer, & ocn_tend_freq_filtered_thickness, & - ocn_tendency_init + ocn_tendency_init, & + ocn_tend_vvel !-------------------------------------------------------------------- ! @@ -235,6 +244,138 @@ subroutine ocn_tend_thick(tendPool, forcingPool)!{{{ end subroutine ocn_tend_thick!}}} +!*********************************************************************** +! +! routine ocn_tend_vvel +! +!> \brief Computes vertical velocity tendency +!> \author Luke Van Roekel +!> \date January 2021 +!> \details +!> This routine computes the vertical velocity tendency for the ocean +! +!----------------------------------------------------------------------- + +subroutine ocn_tend_vvel(tendPool, statePool, forcingPool, meshPool, & + timeLevelIn, dt)!{{{ + + implicit none + + type (mpas_pool_type), intent(inout) :: tendPool !< Input/Output: Tendency structure + type (mpas_pool_type), intent(in) :: statePool !< Input: State information + type (mpas_pool_type), intent(inout) :: forcingPool !< Input: Forcing information + type (mpas_pool_type), intent(in) :: meshPool !< Input: Mesh information + real (kind=RKIND), intent(in) :: dt + integer, intent(in), optional :: timeLevelIn !< Input: Time level for state fields + + type (mpas_pool_type), pointer :: tracersPool + + real (kind=RKIND), dimension(:,:), pointer :: & + normalVelocity, tangentialVelocity, density, potentialDensity, zMid, pressure, & + layerThickness, nonhydrostaticPressure, & + tend_normalVelocity, circulation, relativeVorticity, viscosity, kineticEnergyCell, & + normalizedRelativeVorticityEdge, normalizedPlanetaryVorticityEdge, & + montgomeryPotential, divergence, vertViscTopOfEdge, & + inSituThermalExpansionCoeff, inSituSalineContractionCoeff, vertVelocityNonhydro, & + tend_vertVelocityNonhydro + + ! Scratch Arrays + ! normalThicknessFlux: Flux of thickness through an edge + ! units: m^{2} s^{-1} + real (kind=RKIND), dimension(:,:), allocatable :: normalThicknessFlux + + integer :: timeLevel + + integer :: err, iEdge, iCell, k + integer, pointer :: nVertLevels, indexTemperature, indexSalinity, nEdges, nCells + integer, dimension(:), pointer :: maxLevelCell + + if (present(timeLevelIn)) then + timeLevel = timeLevelIn + else + timeLevel = 1 + end if + + call mpas_pool_get_dimension(meshPool, 'nEdges', nEdges) + call mpas_pool_get_dimension(meshPool, 'nCells', nCells) + + call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) + + call mpas_pool_get_array(statePool, 'vertVelocityNonhydro', vertVelocityNonhydro, timeLevel) + call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity, timeLevel) + call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, timeLevel) + call mpas_pool_get_array(statePool, 'nonhydrostaticPressure', nonhydrostaticPressure, timeLevel) + + call mpas_pool_get_array(tendPool, 'vertVelocityNonhydroTend', tend_vertVelocityNonhydro) + + allocate(normalThicknessFlux(nVertLevels, nEdges+1)) + + ! + ! transport velocity for the tracer. + ! + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(k) + do iEdge = 1, nEdges + do k = 1, nVertLevels + normalThicknessFlux(k, iEdge) = normalTransportVelocity(k, iEdge) * layerThickEdgeFlux(k, iEdge) + end do + end do + !$omp end do + !$omp end parallel + + ! + ! velocity tendency: start accumulating tendency terms + ! + !$omp parallel + !$omp do schedule(runtime) + do iCell = 1, nCells + tend_vertVelocityNonhydro(:, iCell) = 0.0_RKIND + end do + !$omp end do + !$omp end parallel + + if(config_disable_vvel_all_tend) return + + call mpas_timer_start("ocn_tend_vvel") + + call mpas_reconstruct(meshPool, normalVelocity, & + velocityX, velocityY, velocityZ, & + velocityZonal, velocityMeridional, .true.) + ! + ! velocity tendency: nonlinear Coriolis term and grad of kinetic energy + ! + call ocn_vvel_coriolis_tend(meshPool, velocityZonal, tend_vertVelocityNonhydro, err) + + ! + ! velocity tendency: horizontal and vertical advection + ! + call ocn_vvel_advection_tend(meshPool, vertVelocityNonhydro, normalThicknessFlux, & + layerThickness, vertAleTransportTop, dt, & + tend_vertVelocityNonhydro) + + ! + ! velocity tendency: pressure gradient + ! + call ocn_vvel_pgf_tend(meshPool, layerThickness, nonhydrostaticPressure, & + tend_vertVelocityNonhydro, err) + + ! + ! at present, only del2 mixing is implemented for horizontal mixing of + ! vertical velocity + ! + call ocn_vvel_hmix_del2_tend(layerThickEdgeFlux, vertVelocityNonhydro, & + tend_vertVelocityNonhydro, err) + + call ocn_vvel_hmix_del4_tend(layerThickEdgeFlux, vertVelocityNonhydro, & + tend_vertVelocityNonhydro, err) + + + call mpas_timer_stop("ocn_tend_vvel") + deallocate(normalThicknessFlux) + +end subroutine ocn_tend_vvel!}}} + !*********************************************************************** ! ! routine ocn_tend_vel @@ -303,7 +444,9 @@ subroutine ocn_tend_vel(domain, tendPool, statePool, forcingPool, & real (kind=RKIND), dimension(:,:), pointer, contiguous :: & normalVelocity, &! normal velocity - layerThickness + layerThickness, & + vertVelocityNonhydro, & + nonhydrostaticPressure real (kind=RKIND), dimension(:,:,:), pointer, contiguous :: & activeTracers @@ -331,6 +474,10 @@ subroutine ocn_tend_vel(domain, tendPool, statePool, forcingPool, & call mpas_pool_get_array(statePool, 'ssh', ssh, timeLevel) call mpas_pool_get_array(tracersPool, 'activeTracers', & activeTracers, timeLevel) + call mpas_pool_get_array(statePool, 'vertVelocityNonhydro', & + vertVelocityNonhydro, timeLevel) + call mpas_pool_get_array(statePool, 'nonhydrostaticPressure', & + nonhydrostaticPressure, timeLevel) call mpas_pool_get_dimension(tracersPool, 'index_temperature', & indexTemperature) @@ -351,7 +498,8 @@ subroutine ocn_tend_vel(domain, tendPool, statePool, forcingPool, & !$acc density, normRelVortEdge, montgomeryPotential, pressure, & !$acc thermExpCoeff, salineContractCoeff, tangentialVelocity, & !$acc layerThickEdgeFlux, kineticEnergyCell, sfcFlxAttCoeff, potentialDensity, & - !$acc vertAleTransportTop, vertViscTopOfEdge, wettingVelocityFactor) + !$acc vertAleTransportTop, vertViscTopOfEdge, wettingVelocityFactor, & + !$acc nhPressureCorrection) !$acc enter data & !$acc copyin(tendVel, sfcStress, sfcStressMag, & !$acc ssh, normalVelocity, & @@ -410,12 +558,19 @@ subroutine ocn_tend_vel(domain, tendPool, statePool, forcingPool, & layerThickEdgeFlux, normalVelocity, & kineticEnergyCell, tendVel, err) + ! Add nonhydrostatic Coriolis term if nonhydro model is on + if (config_enable_nonhydrostatic_mode) then + call ocn_vel_coriolis_nonhydro_tend(vertVelocityNonhydro, & + tendVel, err) + end if + ! Add vertical advection term -w du/dz call ocn_vel_vadv_tend(normalVelocity, layerThickEdgeFlux, & vertAleTransportTop, tendVel, err) ! Add pressure gradient - call ocn_vel_pressure_grad_tend(ssh, pressure, surfacePressure, & + call ocn_vel_pressure_grad_tend(ssh, nonhydrostaticPressure, & + pressure, surfacePressure, & montgomeryPotential, zMid, & density, potentialDensity, & indxTemp, indxSalt, activeTracers, & diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_del2.F b/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_del2.F index b84e5d39b469..9a7cf717fdb6 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_del2.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_del2.F @@ -29,6 +29,8 @@ module ocn_tracer_hmix_del2 use ocn_constants use ocn_config + use ocn_mesh + use ocn_diagnostics_variables implicit none private @@ -58,7 +60,6 @@ module ocn_tracer_hmix_del2 logical :: del2On real (kind=RKIND) :: eddyDiff2 - !*********************************************************************** contains diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vel_hadv_coriolis.F b/components/mpas-ocean/src/shared/mpas_ocn_vel_hadv_coriolis.F index 7dc0bac48afd..a88c0c457413 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vel_hadv_coriolis.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vel_hadv_coriolis.F @@ -25,6 +25,7 @@ module ocn_vel_hadv_coriolis use ocn_constants use ocn_mesh use ocn_config + use mpas_constants implicit none private @@ -43,6 +44,7 @@ module ocn_vel_hadv_coriolis !-------------------------------------------------------------------- public :: ocn_vel_hadv_coriolis_tend, & + ocn_vel_coriolis_nonhydro_tend, & ocn_vel_hadv_coriolis_init !-------------------------------------------------------------------- @@ -52,8 +54,9 @@ module ocn_vel_hadv_coriolis !-------------------------------------------------------------------- logical :: & - hadvCoriolisDisabled ! on/off switch for hadv/Coriolis tend - ! use disabled since default is on + hadvCoriolisDisabled, & ! on/off switch for hadv/Coriolis tend + ! use disabled since default is on + coriolisNonhydroDisabled logical :: & usePlanetVorticity ! multiplicative mask for including @@ -261,6 +264,68 @@ subroutine ocn_vel_hadv_coriolis_tend(normRelVortEdge, & end subroutine ocn_vel_hadv_coriolis_tend!}}} + subroutine ocn_vel_coriolis_nonhydro_tend(vertVelocityNonhydro, & + tend, err)!{{{ + + !----------------------------------------------------------------- + ! input variables + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(in) :: & + vertVelocityNonhydro !< [in] nonhydro vertical velocity + + !----------------------------------------------------------------- + ! input/output variables + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + tend !< [inout] accumulated velocity tendency + + !----------------------------------------------------------------- + ! output variables + !----------------------------------------------------------------- + + integer, intent(out) :: err !< [out] error flag + + !----------------------------------------------------------------- + ! local variables + !----------------------------------------------------------------- + + integer :: & + iEdge, j, k, &! loop indices for edges, vertical + cell1, cell2, &! neighbor cell indices + kmin, kmax ! min and max indices for vertical + + real (kind=RKIND) :: & + b, & + n1 + + ! End preamble + !----------------------------------------------------------------- + ! Begin code + + !*** Set error code and exit early if disabled + !*** Otherwise, start timer + + err = 0 + if (coriolisNonhydroDisabled) return + + do iEdge = 1, nEdgesOwned + kmin = minLevelEdgeBot(iEdge) + kmax = maxLevelEdgeTop(iEdge) + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + b = 2.0_RKIND * omega * cos(latEdge(iEdge)) + n1 = cos(angleEdge(iEdge)) + do k = kmin, kmax + tend(k,iEdge) = tend(k,iEdge) - & + b*n1*(0.5_RKIND*(vertVelocityNonhydro(k,cell2) & + + vertVelocityNonhydro(k,cell1))) + end do + end do + + end subroutine ocn_vel_coriolis_nonhydro_tend!}}} + !*********************************************************************** ! ! routine ocn_vel_hadv_coriolis_init @@ -291,11 +356,15 @@ subroutine ocn_vel_hadv_coriolis_init(err)!{{{ err = 0 hadvCoriolisDisabled = .false. + coriolisNonhydroDisabled = .false. usePlanetVorticity = .false. !*** Reset module variables based in input configuration - if ( config_disable_vel_coriolis ) hadvCoriolisDisabled = .true. + if ( config_disable_vel_coriolis ) then + hadvCoriolisDisabled = .true. + coriolisNonhydroDisabled = .true. + end if select case (trim(config_time_integrator)) case ('RK4','rk4') diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vel_hmix_del2.F b/components/mpas-ocean/src/shared/mpas_ocn_vel_hmix_del2.F index 4395d6f3d423..e5d2e9398331 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vel_hmix_del2.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vel_hmix_del2.F @@ -24,6 +24,7 @@ module ocn_vel_hmix_del2 use ocn_constants use ocn_config use ocn_mesh + use ocn_diagnostics_variables implicit none private diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vel_pressure_grad.F b/components/mpas-ocean/src/shared/mpas_ocn_vel_pressure_grad.F index 68b547e23202..49008d767ead 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vel_pressure_grad.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vel_pressure_grad.F @@ -93,7 +93,8 @@ module ocn_vel_pressure_grad ! !----------------------------------------------------------------------- - subroutine ocn_vel_pressure_grad_tend(ssh, pressure, surfacePressure, & + subroutine ocn_vel_pressure_grad_tend(ssh, nonhydrostaticPressure, & + pressure, surfacePressure, & montgomeryPotential, zMid, & density, potentialDensity, & indxT, indxS, tracers, & @@ -121,6 +122,9 @@ subroutine ocn_vel_pressure_grad_tend(ssh, pressure, surfacePressure, & density, &!< [in] density potentialDensity !< [in] potentialDensity + real (kind=RKIND), dimension(:,:), intent(in), optional :: & + nonhydrostaticPressure !< Input: nonhydrostatic pressure correction + real (kind=RKIND), dimension(:,:), intent(in) :: & thermExpCoeff, &!< [in] in situ thermal expansion coeff salineContractCoeff !< [in] in situ saline contraction coeff @@ -146,7 +150,7 @@ subroutine ocn_vel_pressure_grad_tend(ssh, pressure, surfacePressure, & !----------------------------------------------------------------- integer :: & - iEdge, k, &! loop indices for edge, cell, vertical loops + iEdge, iCell, k, &! loop indices for edge, cell, vertical loops cell1, cell2, &! neighbor cell indices across edge kMin, kMax ! shallowest and deepest active layer @@ -425,8 +429,8 @@ subroutine ocn_vel_pressure_grad_tend(ssh, pressure, surfacePressure, & !$omp end parallel #endif - !$acc exit data delete(JacobianDxDs) - deallocate(JacobianDxDs) + deallocate(JacobianDxDs) + !$acc exit data delete(JacobianDxDs) case (pGradTypeJacobTS) @@ -586,6 +590,23 @@ subroutine ocn_vel_pressure_grad_tend(ssh, pressure, surfacePressure, & end select ! pGradType + if(config_enable_nonhydrostatic_mode) then + !$omp parallel + !$omp do schedule(runtime) private(cell1, cell2, invdcEdge, k) + do iEdge=1,nEdgesOwned + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + invdcEdge = 1.0_RKIND / dcEdge(iEdge) + + do k=1,maxLevelEdgeTop(iEdge) + tend(k,iEdge) = tend(k,iEdge) + edgeMask(k,iEdge) * invdcEdge * ( & + - ( nonhydrostaticPressure(k,cell2) - nonhydrostaticPressure(k,cell1) ) ) + end do + end do + !$omp end do + !$omp end parallel + endif + call mpas_timer_stop("pressure grad") !-------------------------------------------------------------------- diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vmix.F b/components/mpas-ocean/src/shared/mpas_ocn_vmix.F index 3edb78048919..b703e8d5ba6a 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vmix.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vmix.F @@ -54,6 +54,7 @@ module ocn_vmix public :: ocn_vmix_coefs, & ocn_vel_vmix_tend_implicit, & + ocn_vertVel_vmix_tend_implicit, & ocn_tracer_vmix_tend_implicit, & ocn_vmix_init, & ocn_vmix_implicit @@ -552,6 +553,140 @@ subroutine ocn_vel_vmix_tend_implicit(dt, kineticEnergyCell, vertViscTopOfEdge, end subroutine ocn_vel_vmix_tend_implicit!}}} +!*********************************************************************** +! + ! routine ocn_vertVel_vmix_tend_implicit + ! + !> \brief Computes tendencies for implicit vertical momentum vertical + !mixing + !> \author Luke Van Roekel, Sara Calandrini + !> \date December 2020 + !> \details + !> This routine computes the tendencies for implicit vertical mixing for + !> vertical momentum using computed coefficients. + ! + !----------------------------------------------------------------------- + + subroutine ocn_vertVel_vmix_tend_implicit(meshPool, dt, vertViscTopOfCell, layerThickness, & !{{{ + vertVelocityNonhydro, err) + + implicit none + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + type (mpas_pool_type), intent(in) :: & + meshPool !< Input: mesh information + + real (kind=RKIND), dimension(:,:), intent(in) :: & + vertViscTopOfCell !< Input: vertical mixing coefficients + + real (kind=RKIND), intent(in) :: & + dt !< Input: time step + + real (kind=RKIND), dimension(:,:), intent(in) :: & + layerThickness !< Input: thickness at cell center + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + vertVelocityNonhydro !< Input: velocity + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + integer :: iCell, i, k, N, Nsurf, nCells + integer, pointer :: nVertLevels + integer, dimension(:), pointer :: nCellsArray + + integer, dimension(:), pointer :: maxLevelCell + real(kind=RKIND) :: viscAv, m, A, thickAv + integer, dimension(:,:), pointer :: cellsOnEdge + + real (kind=RKIND), dimension(:), allocatable :: C, velTemp + real (kind=RKIND), dimension(:), allocatable :: bTemp, rTemp + + err = 0 + + if(.not.velVmixOn) return + + call mpas_pool_get_dimension(meshPool, 'nCellsArray', nCellsArray) + call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) + + call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) + call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) + + nCells = nCellsArray( 1 ) + + allocate(bTemp(nVertLevels+1), rTemp(nVertLevels+1)) + allocate(C(nVertLevels+1),velTemp(nVertLevels+1)) + + !$omp parallel + !$omp do schedule(runtime) private(N, Nsurf, C, bTemp, rTemp, k, thickAv, viscAv, A, m, velTemp) + do iCell = 1, nCells + Nsurf = minLevelCell(iCell) + N = maxLevelCell(iCell)+1 + + ! tridiagonal matrix algorithm + C(Nsurf) = 0.0_RKIND ! wnew(Nsurf) = w(Nsurf) + bTemp(Nsurf) = 1.0_RKIND + rTemp(Nsurf) = vertVelocityNonhydro(Nsurf,iCell) + + ! first pass: set the coefficients + do k = Nsurf+1, N-1 + thickAv = 0.5_RKIND * (layerThickness(k-1,iCell) + layerThickness(k,iCell)) + viscAv = 0.5_RKIND * (vertViscTopOfCell(k-1,iCell) + vertViscTopOfCell(k,iCell)) + A = -dt * viscAv / layerThickness(k-1,iCell) / thickAv + m = A / bTemp(k-1) + viscAv = 0.5_RKIND * (vertViscTopOfCell(k+1,iCell) + vertViscTopOfCell(k,iCell)) + C(k) = -dt * viscAv / layerThickness(k,iCell) / thickAv + bTemp(k) = 1.0_RKIND - A - C(k) - m * C(k-1) + rTemp(k) = vertVelocityNonhydro(k,iCell) - m * rTemp(k-1) + end do + + A = 0.0_RKIND ! wnew(N) = w(N) + m = A / bTemp(N-1) + velTemp(N) = (vertVelocityNonhydro(N,iCell) - m * rTemp(N-1)) / (1.0_RKIND - A - m * C(N-1)) + + ! second pass: back substitution + do k = N-1, Nsurf, -1 + velTemp(k) = (rTemp(k) - C(k) * velTemp(k+1)) / bTemp(k) + end do + + vertVelocityNonhydro(Nsurf:N,iCell) = velTemp(Nsurf:N) + vertVelocityNonhydro(1:Nsurf-1,iCell) = 0.0_RKIND + vertVelocityNonhydro(N+1:nVertLevels,iCell) = 0.0_RKIND + + end do + !$omp end do + !$omp end parallel + + deallocate(C,velTemp) + deallocate(bTemp, rTemp) + + !-------------------------------------------------------------------- + + end subroutine ocn_vertVel_vmix_tend_implicit!}}} + + !*********************************************************************** ! ! routine ocn_vel_vmix_tend_implicit_variable @@ -1167,7 +1302,7 @@ subroutine ocn_vmix_implicit(dt, meshPool, statePool, forcingPool, scratchPool, integer :: iCell, timeLevel, k, cell1, cell2, iEdge, iTracer, num_tracers, nCells, nEdges, N, Nsurf real (kind=RKIND), dimension(:), pointer :: bottomDrag, ssh, landIceFraction real (kind=RKIND), dimension(:), allocatable :: landIceEdgeFraction - real (kind=RKIND), dimension(:,:), pointer :: normalVelocity, layerThickness + real (kind=RKIND), dimension(:,:), pointer :: normalVelocity, layerThickness, vertVelocityNonhydro real (kind=RKIND), dimension(:,:), pointer :: tracerGroupSurfaceFlux real (kind=RKIND), dimension(:,:,:), pointer :: tracersGroup real (kind=RKIND), dimension(:,:,:), allocatable :: nonLocalFluxTend @@ -1234,6 +1369,7 @@ subroutine ocn_vmix_implicit(dt, meshPool, statePool, forcingPool, scratchPool, end if call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity, timeLevel) + call mpas_pool_get_array(statePool, 'vertVelocityNonhydro', vertVelocityNonhydro, timeLevel) call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, timeLevel) call mpas_pool_get_array(statePool, 'ssh', ssh, timeLevel) @@ -1355,6 +1491,13 @@ subroutine ocn_vmix_implicit(dt, meshPool, statePool, forcingPool, scratchPool, !$acc exit data copyout(normalVelocity) #endif call mpas_timer_stop('vmix solve momentum') + ! + ! If nonhydrostatic is enabled compute vertical velocity diffusion + ! + if (config_enable_nonhydrostatic_mode) then + call ocn_vertVel_vmix_tend_implicit(meshPool, dt, vertViscTopOfCell, layerThickness, & + vertVelocityNonhydro, err) + endif ! ! Implicit vertical solve for all tracers @@ -1587,6 +1730,7 @@ subroutine ocn_vmix_init(domain, err)!{{{ !-------------------------------------------------------------------- end subroutine ocn_vmix_init!}}} + !*********************************************************************** end module ocn_vmix diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vvel_advection.F b/components/mpas-ocean/src/shared/mpas_ocn_vvel_advection.F new file mode 100644 index 000000000000..4c25b22d9156 --- /dev/null +++ b/components/mpas-ocean/src/shared/mpas_ocn_vvel_advection.F @@ -0,0 +1,167 @@ +! Copyright (c) 2013, Los Alamos National Security, LLC (LANS) +! and the University Corporation for Atmospheric Research (UCAR). +! +! Unless noted otherwise source code is licensed under the BSD license. +! Additional copyright and license information can be found in the LICENSE file +! distributed with this code, or at http://mpas-dev.github.com/license.html +! +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! ocn_vvel_advection +! +!> \brief MPAS ocean vertical velocity advection driver +!> \authors Luke Van Roekel, Darren Engwirda, Sara Calandrini +!> \date November 2021 +!> \details +!> This module contains initialization and driver routines for computing (non-hydrostatic) +!> vertical velocity advection tendencies. It is based, in part, on the tracer advection +!> code. +!> +! +!------------------------------------------------------------------------------- + +module ocn_vvel_advection + + ! module includes + use mpas_kind_types + use mpas_derived_types + use mpas_pool_routines + use mpas_sort + use mpas_hash + use mpas_timer + + use ocn_vvel_horiz_advection + use ocn_vvel_vert_advection + + use ocn_constants + use ocn_config + + implicit none + private + save + + ! public module method interfaces + public :: ocn_vvel_advection_init, & + ocn_vvel_advection_tend + + ! privat module variables + logical :: vvelAdvOn !< flag to turn on tracer advection + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| + + contains + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! routine ocn_vvel_advection_tend +! +!> \brief MPAS ocean vertical velocity advection tendency +!> \author Luke Van Roekel, Darren Engwirda, Sara Calandrini +!> \date November 2021 +!> \details +!> This routine is the driver routine for computing vvel advection +!> tendencies. It simply calls submodule tendency routines based on choice of +!> algorithm. +! +!------------------------------------------------------------------------------- + + subroutine ocn_vvel_advection_tend(meshPool, vertVelocityNonhydro, normalThicknessFlux, & + layerThickness, vertAleTransportTop, dt, & + tend)!{{{ + + implicit none + + !*** Input/Output parameters + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + tend !< [in,out] vertical velocity tendency to which advection added + + !*** Input parameters + + type (mpas_pool_type), intent(in) :: & + meshPool !< [in] mesh information + real (kind=RKIND), dimension(:,:), intent(in) :: & + vertVelocityNonhydro !< [in] non-hydrostatic vertical velocity + real (kind=RKIND), dimension(:,:), intent(in) :: & + normalThicknessFlux !< [in] thickness weighted horizontal velocity + real (kind=RKIND), dimension(:,:), intent(in) :: & + layerThickness !< [in] thickness field at cell centers + real (kind=RKIND), dimension(:,:), intent(in) :: & + vertAleTransportTop !< [in] velocity through layer top levels + real (kind=RKIND), intent(in) :: & + dt !< [in] time-step + + ! end of preamble + !---------------- + ! begin code + + ! immediate return if tracer advection not selected + if(.not. vvelAdvOn) return + + call mpas_timer_start("Vertical Velocity adv") + + call ocn_vvel_horiz_advection_tend(vertVelocityNonhydro, normalThicknessFlux, & + layerThickness, dt, meshPool, tend) + + call ocn_vvel_vert_advection_tend(vertVelocityNonhydro, layerThickness, dt, & + vertAleTransportTop, tend) + + call mpas_timer_stop("Vertical Velocity adv") + + end subroutine ocn_vvel_advection_tend!}}} + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! routine ocn_vvel_advection_init +! +!> \brief MPAS ocean vertical velocity advection tendency +!> \author Luke Van Roekel, Darren Engwirda, Sara Calandrini +!> \date November 2021 +!> \details +!> This routine is the driver for initializing various (non-hydrostatic) +!> vertical velocity advection choices and variables. +! +!------------------------------------------------------------------------------- + + subroutine ocn_vvel_advection_init(err)!{{{ + + implicit none + + !*** output parameters + + integer, intent(out) :: err !< [out] Error flag + + ! end preamble + !------------- + ! begin code + + err = 0 ! initialize error code to success + + ! set some basic flags for options + vvelAdvOn = .not. config_disable_vvel_adv + + ! set all other options from submodule initialization routines + if (err.eq.0) then + call ocn_vvel_horiz_advection_init(err) + end if + + if (err.eq.0) then + call ocn_vvel_vert_advection_init (err) + end if + + ! if an error is returned from init routines, write an error + ! message and return a non-zero error code + if (err /= 0) then + err = 1 + call mpas_log_write( & + 'Error encountered during vertical velocity advection init', & + MPAS_LOG_ERR, masterOnly=.true.) + endif + + end subroutine ocn_vvel_advection_init!}}} + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| + +end module ocn_vvel_advection + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vvel_coriolis.F b/components/mpas-ocean/src/shared/mpas_ocn_vvel_coriolis.F new file mode 100644 index 000000000000..eaf56c2097e4 --- /dev/null +++ b/components/mpas-ocean/src/shared/mpas_ocn_vvel_coriolis.F @@ -0,0 +1,199 @@ +! Copyright (c) 2013, Los Alamos National Security, LLC (LANS) +! and the University Corporation for Atmospheric Research (UCAR). +! +! Unless noted otherwise source code is licensed under the BSD license. +! Additional copyright and license information can be found in the LICENSE file +! distributed with this code, or at http://mpas-dev.github.com/license.html +! +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! ocn_vel_hadv_coriolis +! +!> \brief MPAS ocean tendency of coriolis for vertical momentum equation +!> \author Luke Van Roekel, Darren Engwirda, Sara Calandrini +!> \date December 2020, updated November 2021 +!> \details +!> This module contains the routine for computing +!> tendencies from the vertical momentum coriolis force. +!> +! +!----------------------------------------------------------------------- + +module ocn_vvel_coriolis + + use mpas_timer + use mpas_derived_types + use mpas_pool_routines + use ocn_constants + use ocn_config + use mpas_constants + + implicit none + private + save + + !-------------------------------------------------------------------- + ! + ! Public parameters + ! + !-------------------------------------------------------------------- + + !-------------------------------------------------------------------- + ! + ! Public member functions + ! + !-------------------------------------------------------------------- + + public :: ocn_vvel_coriolis_tend, & + ocn_vvel_coriolis_init + + !-------------------------------------------------------------------- + ! + ! Private module variables + ! + !-------------------------------------------------------------------- + + logical :: vertCoriolisOn + +!*********************************************************************** + +contains + +!*********************************************************************** +! +! routine ocn_vvel_coriolis_tend +! +!> \brief Computes tendency term for coriolis force for vertical momentum equation +!> \author Luke Van Roekel, Darren Engwirda, Sara Calandrini +!> \date January 2021 +!> \details +!> This routine computes the coriolis tendencies for vertical momentum based +!> on current state. +! +!----------------------------------------------------------------------- + + subroutine ocn_vvel_coriolis_tend(meshPool, velocityZonal, tend, err)!{{{ + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(in) :: & + velocityZonal !< Input: Horizontal velocity + + type (mpas_pool_type), intent(in) :: & + meshPool !< Input: mesh information + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + tend !< Input/Output: velocity tendency + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + integer, dimension(:), pointer :: maxLevelCell + + integer :: k, iCell, nCells + integer, pointer :: nVertLevels + integer, dimension(:), pointer :: nCellsArray + real (kind=RKIND) :: fCellExtended + real (kind=RKIND), dimension(:), pointer :: latCell + + err = 0 + + if ( .not. vertCoriolisOn ) return + + call mpas_timer_start("vertical velocity coriolis") + + call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) + + call mpas_pool_get_array(meshPool, 'latCell', latCell) + + call mpas_pool_get_dimension(meshPool, 'nCellsArray', nCellsArray) + call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) + + nCells = nCellsArray( 2 ) + + !FIXME: zonalvelocity may be a bad choice here as it's not updated during a + !timestep but maybe not the worst assumption. + + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(k, iCell, fCellExtended) + do iCell= 1, nCells + + fCellExtended = 2.0_RKIND * omega * cos(latCell(iCell)) + + tend(1,iCell) = tend(1,iCell) + fCellExtended * velocityZonal(1,iCell) + + do k = 2, maxLevelCell(iCell) + tend(k,iCell) = tend(k,iCell) + 0.5_RKIND * & + fCellExtended * (velocityZonal(k,iCell) + velocityZonal(k-1,iCell)) + end do + + end do + !$omp end do + !$omp end parallel + + call mpas_timer_stop("vertical velocity coriolis") + + !-------------------------------------------------------------------- + + end subroutine ocn_vvel_coriolis_tend!}}} + +!*********************************************************************** +! +! routine ocn_vvel_pgf_init +! +!> \brief Initializes ocean vertical momentum PGF tendencies +!> \author Luke Van Roekel +!> \date January 2021 +! +!----------------------------------------------------------------------- + + subroutine ocn_vvel_coriolis_init(err)!{{{ + + !-------------------------------------------------------------------- + + !----------------------------------------------------------------- + ! + ! Output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + err = 0 + + vertCoriolisOn = .true. + if ( config_disable_vvel_coriolis .or. & + .not. config_enable_nonhydrostatic_mode ) vertCoriolisOn = .false. + + !-------------------------------------------------------------------- + + end subroutine ocn_vvel_coriolis_init!}}} + +!*********************************************************************** + +end module ocn_vvel_coriolis + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! vim: foldmethod=marker diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vvel_hmix_del2.F b/components/mpas-ocean/src/shared/mpas_ocn_vvel_hmix_del2.F new file mode 100644 index 000000000000..71840dc9a8d1 --- /dev/null +++ b/components/mpas-ocean/src/shared/mpas_ocn_vvel_hmix_del2.F @@ -0,0 +1,223 @@ +! Copyright (c) 2013, Los Alamos National Security, LLC (LANS) +! and the University Corporation for Atmospheric Research (UCAR). +! +! Unless noted otherwise source code is licensed under the BSD license. +! Additional copyright and license information can be found in the LICENSE file +! distributed with this code, or at http://mpas-dev.github.com/license.html +! +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! ocn_vvel_hmix_del2 +! +!> \brief MPAS ocean horizontal mixing driver for vertical momentum +!> \author Luke VanRoekel, Sara Calandrini +!> \date July 2022 +!> \details +!> This module contains the main driver routine for computing +!> horizontal mixing tendencies in the vertical momentum. +!> +!> It provides an init and a tend function. Each are described below. +! +!----------------------------------------------------------------------- + +module ocn_vvel_hmix_del2 + + use mpas_timer + use mpas_derived_types + use mpas_pool_routines + use mpas_threading + + use ocn_constants + use ocn_config + use ocn_mesh + use ocn_diagnostics_variables + + implicit none + private + save + + !-------------------------------------------------------------------- + ! + ! Public parameters + ! + !-------------------------------------------------------------------- + + !-------------------------------------------------------------------- + ! + ! Public member functions + ! + !-------------------------------------------------------------------- + + public :: ocn_vvel_hmix_del2_tend, & + ocn_vvel_hmix_del2_init + + !-------------------------------------------------------------------- + ! + ! Private module variables + ! + !-------------------------------------------------------------------- + + logical :: del2On + real (kind=RKIND) :: eddyVisc2 + real (kind=RKIND), dimension(:,:), allocatable :: eddyViscArr + + +!*********************************************************************** + +contains + +!*********************************************************************** +! +! routine vvel_hmix_del2_tend +! +!> \brief Computes Laplacian tendency term for horizontal vertical velocity mixing +!> \author Luke Van Roekel -- Copied from tracer_hmix_del2 +!> \date December 2020 +!> \details +!> This routine computes the horizontal mixing tendency for vertical velocity +!> based on current state using a Laplacian parameterization. +! +!----------------------------------------------------------------------- + + subroutine ocn_vvel_hmix_del2_tend(layerThicknessEdge, vertVelocityNonhydro, tend, err)!{{{ + + implicit none + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(in) :: & + layerThicknessEdge !< Input: thickness at edges + + real (kind=RKIND), dimension(:,:), intent(in) :: & + vertVelocityNonhydro !< Input: tracer quantities + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + tend !< Input/Output: velocity tendency + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + integer :: iCell, iEdge, cell1, cell2 + integer :: i, k, nCells, nEdges + + real (kind=RKIND) :: invAreaCell + real (kind=RKIND) :: vvel_turb_flux, flux, r_tmp + + err = 0 + + if (.not.del2On) return + + call mpas_timer_start("vertical velocity del2") + + nCells = nCellsOwned + nEdges = nEdgesHalo( 1 ) + + !allocate(eddyViscArr(nVertLevels+1, nEdges)) + + ! + ! compute a boundary mask to enforce insulating boundary conditions in the horizontal + ! + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(invAreaCell, i, iEdge, cell1, cell2, r_tmp, k, & + !$omp vvel_turb_flux, flux) + do iCell = 1, nCells + invAreaCell = 1.0_RKIND / areaCell(iCell) + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + + r_tmp = meshScalingDel2(iEdge) * eddyVisc2 * dvEdge(iEdge) / dcEdge(iEdge) + + do k = 2, maxLevelEdgeTop(iEdge) + ! \kappa_2 \nabla \phi on edge + vvel_turb_flux = vertVelocityNonhydro(k, cell2) - vertVelocityNonhydro(k, cell1) + + ! div(h \kappa_2 \nabla \phi) at cell center + !flux = 0.5_RKIND*(layerThicknessEdge(k, iEdge) + layerThicknessEdge(k-1,iEdge)) & + ! * vvel_turb_flux * r_tmp + flux = vvel_turb_flux * r_tmp + + tend(k, iCell) = tend(k, iCell) - edgeSignOnCell(i, iCell) * flux * invAreaCell + end do + end do + end do + !$omp end do + !$omp end parallel + + !deallocate(eddyViscArr) + call mpas_timer_stop("vertical velocity del2") + + !-------------------------------------------------------------------- + + end subroutine ocn_vvel_hmix_del2_tend!}}} + +!*********************************************************************** +! +! routine ocn_vvel_hmix_del2_init +! +!> \brief Initializes ocean vertical velocity horizontal mixing quantities +!> \author Luke VanRoekel +!> \date December 2020 +!> \details +!> This routine initializes a variety of quantities related to +!> Laplacian horizontal velocity mixing in the ocean. +! +!----------------------------------------------------------------------- + + subroutine ocn_vvel_hmix_del2_init(err)!{{{ + + !-------------------------------------------------------------------- + + !----------------------------------------------------------------- + ! + ! call individual init routines for each parameterization + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + err = 0 + + del2On = .false. + + !if ( config_use_vertMom_del2 .or. config_use_two_equation_turbulence_model ) then + if ( config_use_vertMom_del2 ) then + if ( config_vertMom_del2 > 0.0_RKIND ) then + del2On = .true. + eddyVisc2 = config_vertMom_del2 + endif + endif + + !-------------------------------------------------------------------- + + end subroutine ocn_vvel_hmix_del2_init!}}} + +!*********************************************************************** + +end module ocn_vvel_hmix_del2 + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! vim: foldmethod=marker diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vvel_hmix_del4.F b/components/mpas-ocean/src/shared/mpas_ocn_vvel_hmix_del4.F new file mode 100644 index 000000000000..0d717f1267b2 --- /dev/null +++ b/components/mpas-ocean/src/shared/mpas_ocn_vvel_hmix_del4.F @@ -0,0 +1,255 @@ +! Copyright (c) 2013, Los Alamos National Security, LLC (LANS) +! and the University Corporation for Atmospheric Research (UCAR). +! +! Unless noted otherwise source code is licensed under the BSD license. +! Additional copyright and license information can be found in the +! LICENSE file +! distributed with this code, or at +! http://mpas-dev.github.com/license.html +! +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! ocn_vvel_hmix_del4 +! +!> \brief MPAS ocean horizontal mixing driver for vertical momentum +!> \author Darren Engwirda, Luke VanRoekel, Sara Calandrini +!> \date April 2023 +!> \details +!> This module contains the main driver routine for computing +!> horizontal mixing tendencies in the vertical momentum. +!> +!> It provides an init and a tend function. Each are described below. +! +!----------------------------------------------------------------------- + +module ocn_vvel_hmix_del4 + + use mpas_timer + use mpas_derived_types + use mpas_pool_routines + use mpas_threading + + use ocn_constants + use ocn_config + use ocn_mesh + use ocn_diagnostics_variables + + implicit none + private + save + + !-------------------------------------------------------------------- + ! + ! Public parameters + ! + !-------------------------------------------------------------------- + + !-------------------------------------------------------------------- + ! + ! Public member functions + ! + !-------------------------------------------------------------------- + + public :: ocn_vvel_hmix_del4_tend, & + ocn_vvel_hmix_del4_init + + !-------------------------------------------------------------------- + ! + ! Private module variables + ! + !-------------------------------------------------------------------- + + logical :: del4On + real (kind=RKIND) :: eddyVisc4 + + +!*********************************************************************** + +contains + +!*********************************************************************** +! +! routine vvel_hmix_del4_tend +! +!> \brief Computes biharmonic tendency term for horizontal vertical +!velocity mixing +!> \author Darren Engwirda, Sara Calandrini +!> \date April 2023 +!> \details +!> This routine computes the horizontal mixing tendency for vertical +!velocity +!> based on current state using a biharmonic parameterization. +! +!----------------------------------------------------------------------- + + subroutine ocn_vvel_hmix_del4_tend(layerThicknessEdge, vertVelocityNonhydro, tend, err)!{{{ + + implicit none + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(in) :: & + layerThicknessEdge !< Input: thickness at edges + + real (kind=RKIND), dimension(:,:), intent(in) :: & + vertVelocityNonhydro !< Input: tracer quantities + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + tend !< Input/Output: velocity tendency + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + integer :: iCell, iEdge, cell1, cell2 + integer :: i, k, nCells, nEdges + + real (kind=RKIND) :: invAreaCell + real (kind=RKIND) :: vvel_turb_flux, flux, r_tmp + real (kind=RKIND), dimension(:, :), allocatable :: del2 + + err = 0 + + if (.not.del4On) return + + call mpas_timer_start("vertical velocity del4") + + nCells = nCellsAll + nEdges = nEdgesHalo( 1 ) + + allocate(del2(nVertLevels+1, nCells)) + del2(:,:) = 0.0_RKIND + + ! + ! compute del^2(vvel) + ! + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(invAreaCell, i, iEdge, cell1, cell2, r_tmp, k, & + !$omp vvel_turb_flux, flux) + do iCell = 1, nCells + invAreaCell = 1.0_RKIND / areaCell(iCell) + !del2(:, iCell) = 0.0_RKIND + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + + r_tmp = dvEdge(iEdge) / dcEdge(iEdge) + + do k = 2, maxLevelEdgeTop(iEdge) + ! \kappa_2 \nabla \phi on edge + vvel_turb_flux = vertVelocityNonhydro(k, cell2) - vertVelocityNonhydro(k, cell1) + + ! div(\kappa_2 \nabla \phi) at cell center + flux = vvel_turb_flux * r_tmp + + del2(k, iCell) = del2(k, iCell) - edgeSignOnCell(i, iCell) * flux * invAreaCell + end do + end do + end do + !$omp end do + + ! + ! compute del^2(del2) + ! + !$omp do schedule(runtime) & + !$omp private(invAreaCell, i, iEdge, cell1, cell2, r_tmp, k, & + !$omp vvel_turb_flux, flux) + do iCell = 1, nCells + invAreaCell = 1.0_RKIND / areaCell(iCell) + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + + r_tmp = meshScalingDel4(iEdge) * eddyVisc4 * dvEdge(iEdge) / dcEdge(iEdge) + + do k = 2, maxLevelEdgeTop(iEdge) + ! \kappa_2 \nabla \phi on edge + vvel_turb_flux = del2(k, cell2) - del2(k, cell1) + + ! div(\kappa_2 \nabla \phi) at cell center + flux = vvel_turb_flux * r_tmp + + ! NB. +ve for del^4 + tend(k, iCell) = tend(k, iCell) + edgeSignOnCell(i, iCell) * flux * invAreaCell + end do + end do + end do + !$omp end do + !$omp end parallel + + deallocate(del2) + + call mpas_timer_stop("vertical velocity del4") + + !-------------------------------------------------------------------- + + end subroutine ocn_vvel_hmix_del4_tend!}}} + +!*********************************************************************** +! +! routine ocn_vvel_hmix_del4_init +! +!> \brief Initializes ocean vertical velocity horizontal mixing quantities +!> \author Darren Engwirda, Sara Calandrini +!> \date April 2023 +!> \details +!> This routine initializes a variety of quantities related to +!> biharmonic vertical velocity mixing in the ocean. +! +!----------------------------------------------------------------------- + + subroutine ocn_vvel_hmix_del4_init(err)!{{{ + + !-------------------------------------------------------------------- + + !----------------------------------------------------------------- + ! + ! call individual init routines for each parameterization + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + err = 0 + + del4On = .false. + + if ( config_use_vertMom_del4 ) then + if ( config_vertMom_del4 > 0.0_RKIND ) then + del4On = .true. + eddyVisc4 = config_vertMom_del4 + endif + endif + + !-------------------------------------------------------------------- + + end subroutine ocn_vvel_hmix_del4_init!}}} + +!*********************************************************************** + +end module ocn_vvel_hmix_del4 + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vvel_horiz_advection.F b/components/mpas-ocean/src/shared/mpas_ocn_vvel_horiz_advection.F new file mode 100644 index 000000000000..cc452f5f2cf1 --- /dev/null +++ b/components/mpas-ocean/src/shared/mpas_ocn_vvel_horiz_advection.F @@ -0,0 +1,187 @@ +! Copyright (c) 2013, Los Alamos National Security, LLC (LANS) +! and the University Corporation for Atmospheric Research (UCAR). +! +! Unless noted otherwise source code is licensed under the BSD license. +! Additional copyright and license information can be found in the LICENSE file +! distributed with this code, or at http://mpas-dev.github.com/license.html +! +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! ocn_vvel_horiz_advection +! +!> \brief MPAS horizontal advection for (non-hydrostatic) vertical velocity +!> \author Luke Van Roekel, Darren Engwirda, Sara Calandrini +!> \date January 2021, updated Novermber 2021 +!> \details +!> This module contains routines for horizontal advection of vertical velocity +!> via a simple TRiSK-type scheme at vertical levels (ie. half-layers). +! +!------------------------------------------------------------------------------- + +module ocn_vvel_horiz_advection + + ! module includes + use mpas_kind_types + use mpas_derived_types + use mpas_pool_routines + use mpas_io_units + use mpas_threading + use ocn_config + use ocn_mesh + use ocn_tracer_advection_shared + use ocn_diagnostics_variables + + implicit none + private + save + + ! public method interfaces + public :: ocn_vvel_horiz_advection_tend, & + ocn_vvel_horiz_advection_init + + contains + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! routine ocn_vvel_horiz_advection_tend +! +!> \brief MPAS horizontal advection tendency for (non-hydrostatic) vertical velocity +!> \author Luke Van Roekel, Darren Engwirda, Sara Calandrini +!> \date 12/21/20, updated 11/19/2021 +!> \details +!> NOTE this is only active for non-hydrostatic mode +! +!----------------------------------------------------------------------- + + subroutine ocn_vvel_horiz_advection_tend(vertVelocityNonhydro, & + normalThicknessFlux, layerThickness, dt, meshPool, tend)!{{{ + + implicit none + + !*** input parameters + + real (kind=RKIND), dimension(:,:), intent(in) :: & + vertVelocityNonhydro !< [in] non-hydrostatic vertical velocity + real (kind=RKIND), dimension(:,:), intent(in) :: & + normalThicknessFlux !< [in] normal thickness flux uh at edges + real (kind=RKIND), dimension(:,:), intent(in) :: & + layerThickness !< [in] layer thickness at cells + real (kind=RKIND), intent(in) :: dt !< [in] time-step + type (mpas_pool_type), intent(in) :: meshPool !< [in] mesh information + + !*** output parameters + real (kind=RKIND), dimension(:,:), intent(inout) :: & + tend !< [out] vvel tendencies + + integer :: iCell, iEdge, jEdge, k, cell1, cell2 + integer :: nVertLevels, nCells, nEdges + real (kind=RKIND) :: vertVelocityFlux + + real (kind=RKIND), dimension(:), allocatable :: invVolCell + real (kind=RKIND), dimension(:,:), allocatable :: horizFlux + + ! Get dimensions + nVertLevels = size(vertVelocityNonhydro,dim=1) + + ! Initialize pointers + nCells = nCellsHalo( 1 ) !1 + nEdges = nEdgesHalo( 2 ) !2 + allocate(invVolCell(nVertLevels)) + allocate(horizFlux(nVertLevels, nEdges)) + + !$omp parallel + !$omp do schedule(runtime) + do iEdge = 1, nEdges + horizFlux(:, iEdge) = 0.0_RKIND + end do + !$omp end do + !$omp end parallel + + ! compute edge fluxes for div(u*h*w) + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(k, cell1, cell2, vertVelocityFlux) + do iEdge = 1, nEdges + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + do k = 1, maxLevelEdgeTop(iEdge) + if (k .gt. 1) then + vertVelocityFlux = & + 0.5_RKIND * (normalThicknessFlux(k-1,iEdge) + & + normalThicknessFlux(k-0,iEdge)) * & + 0.5_RKIND * (vertVelocityNonhydro(k,cell1) + vertVelocityNonhydro(k,cell2)) + else + vertVelocityFlux = & ! surface + 1.0_RKIND * (normalThicknessFlux(k-0,iEdge)) * & + 0.5_RKIND * (vertVelocityNonhydro(k,cell1) + vertVelocityNonhydro(k,cell2)) + end if + + horizFlux(k,iEdge) = & + edgeMask(k,iEdge) * dvEdge(iEdge) * vertVelocityFlux + end do + end do + !$omp end do + !$omp end parallel + + ! eval. tendencies as div(u*h*w) / (A*h) + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(k, iEdge, jEdge, invVolCell) + do iCell = 1, nCells + invVolCell(:) =-1.0e34_RKIND + invVolCell(1) = 1.0_RKIND / (areaCell(iCell) * layerThickness(1,iCell)) + do k = 2, maxLevelCell(iCell) + invVolCell(k) = 2.0_RKIND / ( & + areaCell(iCell) * (layerThickness(k-1,iCell) + layerThickness(k-0,iCell))) + end do + + do jEdge = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(jEdge,iCell) + do k = 1, maxLevelEdgeTop(iEdge) + tend(k,iCell) = tend(k,iCell) + & + edgeSignOnCell(jEdge,iCell) * invVolCell(k) * horizFlux(k,iEdge) + end do + end do + end do + !$omp end do + !$omp end parallel + + deallocate(horizFlux) + deallocate(invVolCell) + + end subroutine ocn_vvel_horiz_advection_tend!}}} + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! routine ocn_vvel_horiz_advection_init +! +!> \brief MPAS initialize vertical velocity advection tendency. +!> \author Luke Van Roekel, Darren Engwirda, Sara Calandrini +!> \date 12/21/2020, updated 11/19/2021 +!> \details +!> This routine initializes constants and choices for the (non-hydrostatic) +!> vertical velocity advection tendency +! +!------------------------------------------------------------------------------- + + subroutine ocn_vvel_horiz_advection_init(err) !{{{ + + implicit none + + !*** input parameters + + !*** output parameters + + integer, intent(out) :: err !< [out] Error Flag + + err = 0 ! set error code to success + + ! don't have actual choices for vvel advection here, unlike tracers + + end subroutine ocn_vvel_horiz_advection_init!}}} + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| + +end module ocn_vvel_horiz_advection + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vvel_pgrad.F b/components/mpas-ocean/src/shared/mpas_ocn_vvel_pgrad.F new file mode 100644 index 000000000000..4b787135e855 --- /dev/null +++ b/components/mpas-ocean/src/shared/mpas_ocn_vvel_pgrad.F @@ -0,0 +1,191 @@ +! Copyright (c) 2013, Los Alamos National Security, LLC (LANS) +! and the University Corporation for Atmospheric Research (UCAR). +! +! Unless noted otherwise source code is licensed under the BSD license. +! Additional copyright and license information can be found in the LICENSE file +! distributed with this code, or at http://mpas-dev.github.com/license.html +! +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! ocn_vvel_pgrad +! +!> \brief MPAS ocean tendency of pressure gradient force for vertical momentum equation +!> \author Luke Van Roekel, Darren Engwirda, Sara Calandrini +!> \date January 2021, updated November 2021 +!> \details +!> This module contains the routine for computing +!> tendencies from the vertical momentum coriolis force. +!> +! +!----------------------------------------------------------------------- + +module ocn_vvel_pgrad + + use mpas_timer + use mpas_derived_types + use mpas_pool_routines + use ocn_constants + use ocn_config + use ocn_mesh + use ocn_diagnostics_variables + + implicit none + private + save + + !-------------------------------------------------------------------- + ! + ! Public parameters + ! + !-------------------------------------------------------------------- + + !-------------------------------------------------------------------- + ! + ! Public member functions + ! + !-------------------------------------------------------------------- + + public :: ocn_vvel_pgf_tend, & + ocn_vvel_pgf_init + + !-------------------------------------------------------------------- + ! + ! Private module variables + ! + !-------------------------------------------------------------------- + + logical :: vertPGFOn + +!*********************************************************************** + +contains + +!*********************************************************************** +! +! routine ocn_vvel_pgf_tend +! +!> \brief Computes tendency term for pressure gradient in vertical momentum equation +!> \author Luke Van Roekel, Darren Engwirda, Sara Calandrini +!> \date January 2021, updated November 2021 +!> \details + +!> This routine computes the PGF tendencies for the vertical momentum based on +!> current state. +! +!----------------------------------------------------------------------- + + subroutine ocn_vvel_pgf_tend(meshPool, layerThickness, nhPressure, tend, err)!{{{ + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(in) :: & + layerThickness !< Input: layer thickness at cells + + real (kind=RKIND), dimension(:,:), intent(in) :: & + nhPressure !< Input: pressure (nonhydrostatic) + + type (mpas_pool_type), intent(in) :: & + meshPool !< Input: mesh information + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + tend !< Input/Output: velocity tendency + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + integer :: j, k, i, iCell, nCells, kmin + integer, pointer :: nVertLevels + integer, dimension(:), pointer :: maxLevelCell, nCellsArray + + err = 0 + + if ( .not. vertPGFOn ) return + + call mpas_timer_start("vertical velocity PGF") + + call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) + call mpas_pool_get_dimension(meshPool, 'nCellsArray', nCellsArray) + call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) + + nCells = nCellsArray( 2 ) + + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(k) + do iCell= 1, nCells + kmin = minLevelCell(iCell) + tend(kmin,iCell) = tend(kmin,iCell) - & ! nhPressure = 0.0 at surface + 2.0_RKIND * (-nhPressure(kmin,iCell)) / layerThickness(kmin,iCell) + do k = kmin+1, maxLevelCell(iCell) + tend(k,iCell) = tend(k,iCell) - & + 2.0_RKIND * (nhPressure(k-1,iCell) - nhPressure(k,iCell)) & + / (layerThickness(k-1,iCell) + layerThickness(k,iCell)) + end do + end do + !$omp end do + !$omp end parallel + + call mpas_timer_stop("vertical velocity PGF") + + !-------------------------------------------------------------------- + + end subroutine ocn_vvel_pgf_tend!}}} + +!*********************************************************************** +! +! routine ocn_vvel_pgf_init +! +!> \brief Initializes ocean momentum PGF tendencies +!> \author Luke Van Roekel +!> \date January 2021 +! +!----------------------------------------------------------------------- + + subroutine ocn_vvel_pgf_init(err)!{{{ + + !-------------------------------------------------------------------- + + !----------------------------------------------------------------- + ! + ! Output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + err = 0 + + vertPGFOn = .true. + if ( .not. config_enable_nonhydrostatic_mode ) vertPGFOn = .false. + + !-------------------------------------------------------------------- + + end subroutine ocn_vvel_pgf_init!}}} + +!*********************************************************************** + +end module ocn_vvel_pgrad + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! vim: foldmethod=marker diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vvel_vert_advection.F b/components/mpas-ocean/src/shared/mpas_ocn_vvel_vert_advection.F new file mode 100644 index 000000000000..e2b48ab1129a --- /dev/null +++ b/components/mpas-ocean/src/shared/mpas_ocn_vvel_vert_advection.F @@ -0,0 +1,135 @@ +! Copyright (c) 2013, Los Alamos National Security, LLC (LANS) +! and the University Corporation for Atmospheric Research (UCAR). +! +! Unless noted otherwise source code is licensed under the BSD license. +! Additional copyright and license information can be found in the LICENSE file +! distributed with this code, or at http://mpas-dev.github.com/license.html +! +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! ocn_vvel_vert_advection +! +!> \brief MPAS vertical advection for (non-hydrostatic) +!> \author Luke Van Roekel, Darren Engwirda, Sara Calandrini +!> \date January 2021 +!> \details +!> This module contains routines for advection of (non-hydrostatic) vertical +!> velocity. +! +!------------------------------------------------------------------------------- + +module ocn_vvel_vert_advection + + ! module includes + use mpas_kind_types + use mpas_derived_types + use mpas_pool_routines + use mpas_io_units + use mpas_threading + use ocn_mesh + use ocn_tracer_advection_shared + + implicit none + private + save + + ! public method interfaces + public :: ocn_vvel_vert_advection_tend, & + ocn_vvel_vert_advection_init + + contains + +!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! routine ocn_vvel_vert_advection_tend +! +!> \brief MPAS vertical advection tendency for (non-hydrostatic) vertical velocity +!> \author Luke Van Roekel, Darren Engwirda, Sara Calandrini +!> \date 12/22/2020, updated 11/20/2021 +!> \details +!> NOTE this is only active for non-hydrostatic mode +! +!----------------------------------------------------------------------- + + subroutine ocn_vvel_vert_advection_tend(vertVelocityNonhydro, layerThickness, dt, & + vertAleTransportTop, tend)!{{{ + + implicit none + + !*** input parameters + + real (kind=RKIND), dimension(:,:), intent(in) :: & + vertVelocityNonhydro !< Input: non-hydrostatic vertical velocity + real (kind=RKIND), dimension(:,:), intent(in) :: & + layerThickness !< Input: layer thickness at cells + real (kind=RKIND), intent(in) :: dt !< Input: time-step + real (kind=RKIND), dimension(:,:), intent(in) :: & + vertAleTransportTop !< Input: velocity through layer top levels + + !*** output parameters + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + tend !< Output: vvel tendencies + + integer :: iCell, k + real (kind=RKIND) :: levelThickness, vertVelocityUp, vertVelocityDn + + ! -w_ale * dw/dz + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(k, levelThickness, vertVelocityUp, vertVelocityDn) + do iCell = 1, nCellsOwned + ! w_ale = 0.0 at k==1 + do k = 2, maxLevelCell(iCell) + levelThickness = & + 0.5_RKIND * (layerThickness(k-1,iCell) + layerThickness(k,iCell)) + + vertVelocityUp = & + 0.5_RKIND * (vertVelocityNonhydro(k-1,iCell) + vertVelocityNonhydro(k,iCell)) + vertVelocityDn = & + 0.5_RKIND * (vertVelocityNonhydro(k+1,iCell) + vertVelocityNonhydro(k,iCell)) + + tend(k,iCell) = tend(k,iCell) - & + vertAleTransportTop(k,iCell) * & + (vertVelocityUp - vertVelocityDn) / levelThickness + end do + end do + !$omp end do + !$omp end parallel + + end subroutine ocn_vvel_vert_advection_tend!}}} + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! routine ocn_vvel_vert_advection_init +! +!> \brief MPAS initialize vertical velocity advection tendency. +!> \author Luke Van Roekel, Darren Engwirda, Sara Calandrini +!> \date January 2021, updated November 2021 +!> \details +!> This routine initializes constants and choices for the (non-hydrostatic) +!> vertical velocity advection tendency +! +!------------------------------------------------------------------------------- + + subroutine ocn_vvel_vert_advection_init(err) !{{{ + + implicit none + + !*** input parameters + + !*** output parameters + + integer, intent(out) :: err !< [out] Error Flag + + err = 0 ! set error code to success + + ! don't have actual choices for vvel advection here, unlike tracers + + end subroutine ocn_vvel_vert_advection_init!}}} + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| + +end module ocn_vvel_vert_advection + +!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||