diff --git a/components/homme/cmake/machineFiles/greatlakes.cmake b/components/homme/cmake/machineFiles/greatlakes.cmake new file mode 100644 index 000000000000..6328227fd8df --- /dev/null +++ b/components/homme/cmake/machineFiles/greatlakes.cmake @@ -0,0 +1,37 @@ +# CMake initial cache file for Linux 64bit RHEL7/CENTOS7 +# tested with stock gcc/gfortran and packages from EPEL: +# openmpi-devel +# blas-devel +# lapack-devel +# netcdf-fortran-devel +# +# +SET (CMAKE_Fortran_COMPILER mpif90 CACHE FILEPATH "") +SET (CMAKE_C_COMPILER mpicc CACHE FILEPATH "") +SET (CMAKE_CXX_COMPILER mpicc CACHE FILEPATH "") +SET (NetCDF_C $ENV{NETCDF_C_PATH} CACHE FILEPATH "") +SET (NetCDF_C_LIBRARY $ENV{NETCDF_C_PATH}/lib/libnetcdf.so CACHE FILEPATH "") +SET (NetCDF_C_INCLUDE_DIR $ENV{NETCDF_C_PATH}/include CACHE FILEPATH "") +SET (NetCDF_Fortran $ENV{NETCDF_F_PATH} CACHE FILEPATH "") + +SET (NetCDF_Fortran_PATH $ENV{NETCDF_F_PATH} CACHE FILEPATH "") +SET (NetCDF_Fortran_LIBRARY $ENV{NETCDF_F_PATH}/lib/libnetcdff.so CACHE FILEPATH "") +SET (NetCDF_Fortran_INCLUDE_DIR $ENV{NETCDF_F_PATH}/include CACHE FILEPATH "") +SET (HDF5_C_LIBRARY $ENV{HDF5_PATH}/lib/libhdf5.so CACHE FILEPATH "") +SET (HDF5_C_INCLUDE_DIR $ENV{HDF5_PATH}/include CACHE FILEPATH "") +SET (HDF5_HL_LIBRARY $ENV{HDF5_PATH}/lib/libhdf5_hl.so CACHE FILEPATH "") +SET (HDF5_HL_INCLUDE_DIR $ENV{HDF5_PATH}/include CACHE FILEPATH "") +SET (PNETCDF_DIR $ENV{PNETCDF_PATH} CACHE FILEPATH "") + + + + + + +SET (WITH_PNETCDF TRUE CACHE FILEPATH "") + +# hack until findnetcdf is updated to look for netcdf.mod +#SET (ADD_Fortran_FLAGS $ENV{NETCDF_F_PATH}/include CACHE STRING "") + +SET (USE_QUEUING FALSE CACHE BOOL "") +SET (HOMME_FIND_BLASLAPACK ON CACHE BOOL "") diff --git a/components/homme/src/arkode/theta-l/arkode_interface.F90 b/components/homme/src/arkode/theta-l/arkode_interface.F90 index f28d61c0ca14..e2378b0def40 100644 --- a/components/homme/src/arkode/theta-l/arkode_interface.F90 +++ b/components/homme/src/arkode/theta-l/arkode_interface.F90 @@ -248,7 +248,7 @@ function FColumnSolSolve(b, t, y, gamma, x) result(ierr) use kinds, only: real_kind use HommeNVector, only: NVec_t use hybvcoord_mod, only: hvcoord_t - use physical_constants, only: g + use physical_constants, only: gravit use dimensions_mod, only: np, nlev, nlevp use parallel_mod, only: abortmp use iso_c_binding diff --git a/components/homme/src/implicit_mod.F90 b/components/homme/src/implicit_mod.F90 index dce285c230aa..af7d87ad44b7 100644 --- a/components/homme/src/implicit_mod.F90 +++ b/components/homme/src/implicit_mod.F90 @@ -9,7 +9,7 @@ subroutine advance_imp_nonstag(elem, edge1, edge2, edge3, red, deriv, cg, & use kinds, only : real_kind use dimensions_mod, only : np, nlev, nvar - use physical_constants, only : g + use physical_constants, only : gravit use element_mod, only : element_t use parallel_mod, only : parallel_t, abortmp use edge_mod, only : edgevpack, edgevunpack @@ -530,7 +530,7 @@ subroutine precon_gmres(vv, z, nelemd, xstate, c_ptr_to_object, c_ptr_to_pre) & use ,intrinsic :: iso_c_binding use kinds, only : real_kind use dimensions_mod, only : np, nlev, nvar, nelem - use physical_constants, only : g + use physical_constants, only : gravit use element_mod, only : element_t ! use parallel_mod, only : parallel_t use edge_mod, only : edgevpack, edgevunpack @@ -630,7 +630,7 @@ subroutine precon_si(vv, z, nelemd, xstate, c_ptr_to_object, c_ptr_to_pre) & use ,intrinsic :: iso_c_binding use kinds, only : real_kind - use physical_constants, only : rearth, g + use physical_constants, only : rearth, gravit use dimensions_mod, only : np, nlev, nvar, nelem use element_mod, only : element_t use edge_mod, only : edgevpack, edgevunpack @@ -920,7 +920,7 @@ subroutine sw_picard_schur_op(xs, nelemd,fx, c_ptr_to_object) bind(C,name='sw_pi use bndry_mod, only : bndry_exchangev use derived_type_mod, only : derived_type use perf_mod, only : t_startf, t_stopf - use physical_constants, only : g + use physical_constants, only : gravit implicit none real (c_double) ,intent(in) :: xs(nelemd) @@ -1227,7 +1227,7 @@ subroutine sw_picard_DFinvBt_op(xs, nelemd,fx,nret, c_ptr_to_object) bind(C,name use bndry_mod, only : bndry_exchangev use derived_type_mod, only : derived_type use perf_mod, only : t_startf, t_stopf - use physical_constants, only : g + use physical_constants, only : gravit implicit none real (c_double) ,intent(in) :: xs(nelemd) @@ -1405,7 +1405,7 @@ subroutine sw_picard_block_11(xs, nelemd,fx, c_ptr_to_object) bind(C,name='sw_pi use bndry_mod, only : bndry_exchangev use derived_type_mod, only : derived_type use perf_mod, only : t_startf, t_stopf - use physical_constants, only : g + use physical_constants, only : gravit implicit none real (c_double) ,intent(in) :: xs(nelemd) @@ -1615,7 +1615,7 @@ subroutine sw_picard_block_21(xs, nelemd,fx, nret, c_ptr_to_object) bind(C,name= use bndry_mod, only : bndry_exchangev use derived_type_mod, only : derived_type use perf_mod, only : t_startf, t_stopf - use physical_constants, only : g + use physical_constants, only : gravit implicit none real (c_double) ,intent(in) :: xs(nelemd) @@ -1789,7 +1789,7 @@ subroutine sw_jacobian_op(xs, nelemd,fx, c_ptr_to_object) bind(C,name='sw_jacobi use bndry_mod, only : bndry_exchangev use derived_type_mod, only : derived_type use perf_mod, only : t_startf, t_stopf - use physical_constants, only : g + use physical_constants, only : gravit implicit none real (c_double) ,intent(in) :: xs(nelemd) diff --git a/components/homme/src/interp_movie_mod.F90 b/components/homme/src/interp_movie_mod.F90 index e4e1bf258035..8a9230280c41 100644 --- a/components/homme/src/interp_movie_mod.F90 +++ b/components/homme/src/interp_movie_mod.F90 @@ -30,7 +30,7 @@ module interp_movie_mod nf_init_decomp, & get_varindex - use physical_constants, only : omega, g, rrearth, dd_pi, kappa, p0 + use physical_constants, only : omega, gravit, rrearth, dd_pi, kappa, p0 use control_mod, only : test_case, runtype, & restartfreq, & diff --git a/components/homme/src/prim_movie_mod.F90 b/components/homme/src/prim_movie_mod.F90 index 34279bbfb1dd..79179e363721 100644 --- a/components/homme/src/prim_movie_mod.F90 +++ b/components/homme/src/prim_movie_mod.F90 @@ -50,7 +50,7 @@ module prim_movie_mod nf_get_frame use coordinate_systems_mod, only : cartesian2D_t, spherical_polar_t, cartesian3D_t, spherical_to_cart - use physical_constants, only : g, kappa, p0, dd_pi + use physical_constants, only : gravit, kappa, p0, dd_pi use dof_mod, only : UniquePoints, UniqueCoords, UniqueNcolsP, createmetadata #ifndef HOMME_WITHOUT_PIOLIBRARY use pio, only : io_desc_t, pio_iotask_rank !_EXTERNAL diff --git a/components/homme/src/share/control_mod.F90 b/components/homme/src/share/control_mod.F90 index c402fc7cf956..c6b1455f6f01 100644 --- a/components/homme/src/share/control_mod.F90 +++ b/components/homme/src/share/control_mod.F90 @@ -219,6 +219,8 @@ module control_mod real (kind=real_kind), public :: pertlim = 0 !pertibation to temperature [like CESM] #endif + logical, public :: atm_is_deep = .true. ! height of mountain range (meters) + logical, public :: gravity_is_const = .false. ! height of mountain range (meters) ! shallow water advection test paramters ! kmass = level index with density. other levels contain test tracers integer, public :: kmass = -1 diff --git a/components/homme/src/share/cube_mod.F90 b/components/homme/src/share/cube_mod.F90 index 2d8fb1293d25..13dd2b47767a 100644 --- a/components/homme/src/share/cube_mod.F90 +++ b/components/homme/src/share/cube_mod.F90 @@ -774,6 +774,7 @@ end subroutine dmap_elementlocal subroutine coriolis_init_atomic(elem) use element_mod, only : element_t use physical_constants, only : omega + use control_mod, only: atm_is_deep type (element_t) :: elem @@ -785,7 +786,10 @@ subroutine coriolis_init_atomic(elem) rangle = rotate_grid*DD_PI/180 do j=1,np do i=1,np - if ( rotate_grid /= 0) then + if ( rotate_grid /= 0 ) then + if (atm_is_deep) then + call abortmp('Rotated grid not yet supported in deep atmosphere') + end if lat = elem%spherep(i,j)%lat lon = elem%spherep(i,j)%lon elem%fcor(i,j)= 2*omega* & diff --git a/components/homme/src/share/deep_atm_ctrl_mod.F90 b/components/homme/src/share/deep_atm_ctrl_mod.F90 new file mode 100644 index 000000000000..640721179c5f --- /dev/null +++ b/components/homme/src/share/deep_atm_ctrl_mod.F90 @@ -0,0 +1,248 @@ +module deep_atm_mod +use physical_constants, only : rearth, gravit, omega +use kinds, only: real_kind +use dimensions_mod, only: np +implicit none +private +save + +public :: r_hat_from_phi, & + z_from_phi, & + g_from_phi, & + phi_from_z, & + quasi_hydrostatic_terms +public :: atm_is_deep, gravity_is_const + +logical, parameter :: atm_is_deep=.true. +logical, parameter :: gravity_is_const=.false. +interface r_hat_from_phi + module procedure r_hat_from_phi_rank_ijk + module procedure r_hat_from_phi_rank_ij + module procedure r_hat_from_phi_rank_k + module procedure r_hat_from_phi_no_rank +end interface +interface z_from_phi + module procedure z_from_phi_rank_ijk + module procedure z_from_phi_rank_ij + module procedure z_from_phi_rank_k + module procedure z_from_phi_no_rank +end interface +interface g_from_phi + module procedure g_from_phi_rank_ijk + module procedure g_from_phi_rank_ij + module procedure g_from_phi_rank_k + module procedure g_from_phi_no_rank +end interface +interface phi_from_z + module procedure phi_from_z_rank_ijk + module procedure phi_from_z_rank_ij + module procedure phi_from_z_rank_k + module procedure phi_from_z_no_rank +end interface + + + + +contains +function r_hat_from_phi_rank_ijk(phi, k) result(r_hat) + integer, intent(in) :: k + real(kind=real_kind), intent(in) :: phi(np,np,k) + real(kind=real_kind) :: r_hat(np,np,k) + r_hat = elementwise_ijk(phi, k, "r_hat_from_phi_scalar") +end function r_hat_from_phi_rank_ijk +function r_hat_from_phi_rank_ij(phi) result(r_hat) + real(kind=real_kind), intent(in) :: phi(np,np) + real(kind=real_kind) :: r_hat(np,np) + r_hat = elementwise_ij(phi, "r_hat_from_phi_scalar") +end function r_hat_from_phi_rank_ij +function r_hat_from_phi_rank_k(phi, k) result(r_hat) + integer, intent(in) :: k + real(kind=real_kind), intent(in) :: phi(k) + real(kind=real_kind) :: r_hat(k) + r_hat = elementwise_k(phi, k, "r_hat_from_phi_scalar") +end function r_hat_from_phi_rank_k +function r_hat_from_phi_no_rank(phi) result(r_hat) + real(kind=real_kind), intent(in) :: phi + real(kind=real_kind) :: r_hat + r_hat = r_hat_from_phi_scalar(phi) +end function r_hat_from_phi_no_rank +function z_from_phi_rank_ijk(phi, k) result(z) + integer, intent(in) :: k + real(kind=real_kind), intent(in) :: phi(np,np,k) + real(kind=real_kind) :: z(np,np,k) + z = elementwise_ijk(phi, k, "z_from_phi_scalar") +end function z_from_phi_rank_ijk +function z_from_phi_rank_ij(phi) result(z) + real(kind=real_kind), intent(in) :: phi(np,np) + real(kind=real_kind) :: z(np,np) + z = elementwise_ij(phi, "z_from_phi_scalar") +end function z_from_phi_rank_ij +function z_from_phi_rank_k(phi, k) result(z) + integer, intent(in) :: k + real(kind=real_kind), intent(in) :: phi(k) + real(kind=real_kind) :: z(k) + z = elementwise_k(phi, k, "z_from_phi_scalar") +end function z_from_phi_rank_k +function z_from_phi_no_rank(phi) result(z) + real(kind=real_kind), intent(in) :: phi + real(kind=real_kind) :: z + z = z_from_phi_scalar(phi) +end function z_from_phi_no_rank +function phi_from_z_rank_ijk(z, k) result(phi) + integer, intent(in) :: k + real(kind=real_kind), intent(in) :: z(np,np,k) + real(kind=real_kind) :: phi(np,np,k) + phi = elementwise_ijk(z, k, "phi_from_z_scalar") +end function phi_from_z_rank_ijk +function phi_from_z_rank_ij(z) result(phi) + real(kind=real_kind), intent(in) :: z(np,np) + real(kind=real_kind) :: phi(np,np) + phi = elementwise_ij(z, "phi_from_z_scalar") +end function phi_from_z_rank_ij +function phi_from_z_rank_k(z, k) result(phi) + integer, intent(in) :: k + real(kind=real_kind), intent(in) :: z(k) + real(kind=real_kind) :: phi(k) + phi = elementwise_k(z, k, "phi_from_z_scalar") +end function phi_from_z_rank_k +function phi_from_z_no_rank(z) result(phi) + real(kind=real_kind), intent(in) :: z + real(kind=real_kind) :: phi + phi = phi_from_z_scalar(z) +end function phi_from_z_no_rank +function g_from_phi_rank_ijk(phi, k) result(g) + integer, intent(in) :: k + real(kind=real_kind), intent(in) :: phi(np,np,k) + real(kind=real_kind) :: g(np,np,k) + g = elementwise_ijk(phi, k, "g_from_phi_scalar") +end function g_from_phi_rank_ijk +function g_from_phi_rank_ij(phi) result(g) + real(kind=real_kind), intent(in) :: phi(np,np) + real(kind=real_kind) :: g(np,np) + g = elementwise_ij(phi, "g_from_phi_scalar") +end function g_from_phi_rank_ij +function g_from_phi_rank_k(phi, k) result(g) + integer, intent(in) :: k + real(kind=real_kind), intent(in) :: phi(k) + real(kind=real_kind) :: g(k) + g = elementwise_k(phi, k, "g_from_phi_scalar") +end function g_from_phi_rank_k +function g_from_phi_no_rank(phi) result(g) + real(kind=real_kind), intent(in) :: phi + real(kind=real_kind) :: g + g = g_from_phi_scalar(phi) +end function g_from_phi_no_rank +function dispatch(inval, funcname) result(outval) + real(kind=real_kind), intent(in) :: inval + character(*), intent(in) :: funcname + real(kind=real_kind) :: outval + + select case (funcname) + case ("r_hat_from_phi_scalar") + outval = r_hat_from_phi_scalar(inval) + case ("z_from_phi_scalar") + outval = z_from_phi_scalar(inval) + case ("g_from_phi_scalar") + outval = g_from_phi_scalar(inval) + case ("phi_from_z_scalar") + outval = phi_from_z_scalar(inval) + end select +end function dispatch + +function elementwise_ijk(inval, k, funcname) result(outval) + integer, intent(in) :: k + real(kind=real_kind), intent(in) :: inval(np,np,k) + real(kind=real_kind) :: outval(np,np,k) + character (*), intent(in) :: funcname + + integer :: i, j, kk + + do kk=1,size(inval, 3) + do j=1,size(inval, 2) + do i=1,size(inval,1) + outval(i,j,kk) = dispatch(inval(i, j, kk), funcname) + end do + end do + end do +end function elementwise_ijk +function elementwise_ij(inval, funcname) result(outval) + real(kind=real_kind), intent(in) :: inval(np,np) + real(kind=real_kind) :: outval(np,np) + character (*), intent(in) :: funcname + integer :: i, j + + do j=1,size(inval, 2) + do i=1,size(inval,1) + outval(i,j) = dispatch(inval(i, j), funcname) + end do + end do +end function elementwise_ij +function elementwise_k(inval, k, funcname) result(outval) + integer, intent(in) :: k + real(kind=real_kind), intent(in) :: inval(k) + real(kind=real_kind) :: outval(k) + character (*), intent(in) :: funcname + + integer :: kk + + do kk=1,size(inval, 1) + outval(kk) = dispatch(inval(kk), funcname) + end do +end function elementwise_k + + +function g_from_z_scalar(z) result(g) + real(kind=real_kind), intent(in) :: z + real(kind=real_kind) :: g + if (atm_is_deep .and. .not. gravity_is_const) then + g = gravit / ((rearth + z)/rearth)**2 + else + g = gravit + end if +end function g_from_z_scalar + +function r_hat_from_phi_scalar(phi) result(r_hat) + real(kind=real_kind), intent(in) :: phi + real(kind=real_kind) :: r_hat + if (atm_is_deep) then + r_hat = (rearth + z_from_phi_scalar(phi))/rearth + else + r_hat = 1.0_real_kind + end if + +end function r_hat_from_phi_scalar + +function z_from_phi_scalar(phi) result(z) + real(kind=real_kind), intent(in) :: phi + real(kind=real_kind) :: z + real(kind=real_kind) :: b + if (atm_is_deep .and. .not. gravity_is_const) then + b = (2 * phi * rearth -gravit*rearth**2) + z = -2*phi*rearth**2/(b-sqrt(b**2 - 4 * phi**2 * rearth**2) ) + else + z = phi/gravit + end if +end function z_from_phi_scalar + +function phi_from_z_scalar(z) result(phi) + real(kind=real_kind), intent(in) :: z + real(kind=real_kind) :: phi + phi = z * g_from_z_scalar(z) +end function phi_from_z_scalar + +function g_from_phi_scalar(phi) result(g) + real(kind=real_kind), intent(in) :: phi + real(kind=real_kind) :: g + g = g_from_z_scalar(z_from_phi_scalar(phi)) +end function g_from_phi_scalar + + +function quasi_hydrostatic_terms(u_top, lat, phi_top) result(terms) + real(kind=real_kind), intent(in) :: u_top(np,np,2) + real(kind=real_kind), intent(in) :: lat(np,np) + real(kind=real_kind), intent(in) :: phi_top(np,np) + real(kind=real_kind) :: terms(np,np) + terms = - (u_top(:,:,1)**2 + u_top(:,:,2)**2)/(rearth + z_from_phi(phi_top)) - 2 * omega * cos(lat) * u_top(:,:,1) + +end function quasi_hydrostatic_terms +end module diff --git a/components/homme/src/share/gllfvremap_mod.F90 b/components/homme/src/share/gllfvremap_mod.F90 index 1e0ee9b81842..2d82314af3cd 100644 --- a/components/homme/src/share/gllfvremap_mod.F90 +++ b/components/homme/src/share/gllfvremap_mod.F90 @@ -556,7 +556,8 @@ end subroutine gfr_dyn_to_fv_phys_topo_data subroutine gfr_dyn_to_fv_phys_topo_data_elem(ie, elem, square, augment_variance, g, p) ! Element-level impl of gfr_dyn_to_fv_phys_topo_data. - use physical_constants, only: grav => g + use physical_constants, only: gravit + use deep_atm_mod, only: z_from_phi integer, intent(in) :: ie type (element_t), intent(in) :: elem(:) @@ -578,7 +579,7 @@ subroutine gfr_dyn_to_fv_phys_topo_data_elem(ie, elem, square, augment_variance, do k = 1,nf2 ! Integrate (phis_gll - phis_fv)^2 over FV subcell (i,j). Do this ! using gfr_g2f_scalar; thus, only one entry out of nf^2 is used. - wg(:,:,1) = ((elem(ie)%state%phis - phispg(k))/grav)**2 + wg(:,:,1) = (z_from_phi(elem(ie)%state%phis) - z_from_phi(phispg(k)))**2 ! DA_CHANGE call gfr_g2f_scalar(ie, elem(ie)%metdet, wg(:,:,:1), wf(:,1:1)) ! Use just entry (i,j). wf(k,2) = max(zero, wf(k,1)) diff --git a/components/homme/src/share/gllfvremap_util_mod.F90 b/components/homme/src/share/gllfvremap_util_mod.F90 index 4e56f1a514a4..86742e121a00 100644 --- a/components/homme/src/share/gllfvremap_util_mod.F90 +++ b/components/homme/src/share/gllfvremap_util_mod.F90 @@ -82,7 +82,7 @@ end subroutine finish subroutine set_state(s, nt1, nt2, ntq, elem) ! Convenience wrapper to set_elem_state. - use physical_constants, only: g + use physical_constants, only: gravit use element_ops, only: set_elem_state type (State_t), intent(in) :: s @@ -90,7 +90,7 @@ subroutine set_state(s, nt1, nt2, ntq, elem) type (element_t), intent(inout) :: elem call set_elem_state(s%u, s%v, s%w, s%wi, s%T, s%ps, s%phis, s%p, s%dp, s%z, s%zi, & - g, elem, nt1, nt2, ntq) + gravit, elem, nt1, nt2, ntq) !DA_CHANGE end subroutine set_state subroutine set_gll_state(hvcoord, elem, nt1, nt2) @@ -102,10 +102,11 @@ subroutine set_gll_state(hvcoord, elem, nt1, nt2) ! function on the sphere. use dimensions_mod, only: nlev, qsize - use physical_constants, only: g, dd_pi + use physical_constants, only: gravit, dd_pi use coordinate_systems_mod, only: cartesian3D_t, change_coordinates use hybvcoord_mod, only: hvcoord_t use element_ops, only: get_field + use deep_atm_mod, only: g_from_phi type (hvcoord_t) , intent(in) :: hvcoord type (element_t), intent(inout) :: elem @@ -156,7 +157,7 @@ subroutine set_gll_state(hvcoord, elem, nt1, nt2) ! a bit of a kludge call set_state(s1, nt1, nt2, nt1, elem) call get_field(elem, 'rho', wr(:,:,:,1), hvcoord, nt1, nt1) - s1%w = -elem%derived%omega_p/(wr(:,:,:,1)*g) + s1%w = -elem%derived%omega_p/(wr(:,:,:,1)*gravit) ! WARNING: NOT ADAPTED TO DEEP ATMOSPHERE s1%wi(:,:,:nlev) = s1%w s1%wi(:,:,nlevp) = s1%w(:,:,nlev) call set_state(s1, nt1, nt2, nt1, elem) @@ -182,7 +183,7 @@ function run(hybrid, hvcoord, elem, nets, nete, nphys, tendency) result(nerr) use parallel_mod, only: global_shared_buf, global_shared_sum use global_norms_mod, only: wrap_repro_sum use reduction_mod, only: ParallelMin, ParallelMax - use physical_constants, only: g, p0, kappa + use physical_constants, only: gravit, p0, kappa use edge_mod, only: edgevpack_nlyr, edgevunpack_nlyr, edge_g use bndry_mod, only: bndry_exchangev use control_mod, only: ftype diff --git a/components/homme/src/share/namelist_mod.F90 b/components/homme/src/share/namelist_mod.F90 index dfe761e206a0..717fb264a61a 100644 --- a/components/homme/src/share/namelist_mod.F90 +++ b/components/homme/src/share/namelist_mod.F90 @@ -57,6 +57,8 @@ module namelist_mod prescribed_wind, & ftype, & limiter_option,& + atm_is_deep, & + gravity_is_const, & nu, & nu_s, & nu_q, & @@ -282,6 +284,8 @@ subroutine readnl(par) disable_diagnostics, & prescribed_wind, & se_ftype, & ! forcing type + atm_is_deep, & + gravity_is_const, & nu, & nu_s, & nu_q, & @@ -788,6 +792,10 @@ subroutine readnl(par) call MPI_bcast(hv_theta_thresh,1, MPIreal_t, par%root,par%comm,ierr) call MPI_bcast(vert_remap_q_alg,1, MPIinteger_t, par%root,par%comm,ierr) call MPI_bcast(vert_remap_u_alg,1, MPIinteger_t, par%root,par%comm,ierr) + + + call MPI_bcast(gravity_is_const, 1, MPIlogical_t , par%root,par%comm,ierr) + call MPI_bcast(atm_is_deep, 1, MPIlogical_t , par%root,par%comm,ierr) call MPI_bcast(nu, 1, MPIreal_t , par%root,par%comm,ierr) call MPI_bcast(nu_s, 1, MPIreal_t , par%root,par%comm,ierr) diff --git a/components/homme/src/share/physical_constants.F90 b/components/homme/src/share/physical_constants.F90 index b4c9aa10c47a..17c4ce9f6ff0 100644 --- a/components/homme/src/share/physical_constants.F90 +++ b/components/homme/src/share/physical_constants.F90 @@ -13,7 +13,7 @@ module physical_constants ! set physical constants from CAM variables use physconst, only: & pi, & ! _EXTERNAL - g => gravit, & + gravit=>gravit, & rearth, & omega, & Rgas => rair, & @@ -38,7 +38,7 @@ module physical_constants real (kind=real_kind), public, parameter :: DD_PI = pi real (kind=longdouble_kind), public, parameter :: QQ_PI = 3.141592653589793238462643383279_longdouble_kind public :: rearth ! m - public :: g ! m s^-2 + public :: gravit ! m s^-2 public :: omega ! s^-1 public :: Rgas real (kind=real_kind), public, parameter :: Cp = cpair @@ -57,8 +57,8 @@ module physical_constants real (kind=longdouble_kind), public, parameter :: QQ_PI = 3.141592653589793238462643383279_longdouble_kind real (kind=real_kind), public, parameter :: rearth0 = 6.376D6 ! m real (kind=real_kind), public :: rearth = rearth0 ! m - real (kind=real_kind), public, parameter :: g = 9.80616D0 ! m s^-2 - real (kind=real_kind), public :: ginv = 1.0_real_kind/g + real (kind=real_kind), public, parameter :: gravit = 9.80616D0 ! m s^-2 + real (kind=real_kind), public :: gravitinv = 1.0_real_kind/gravit real (kind=real_kind), public, parameter :: omega0 = 7.292D-5 ! s^-1 real (kind=real_kind), public :: omega = omega0 real (kind=real_kind), public, parameter :: Rgas = 287.04D0 diff --git a/components/homme/src/share/physics_mod.F90 b/components/homme/src/share/physics_mod.F90 index e954360c8f36..8dc623bcfa55 100644 --- a/components/homme/src/share/physics_mod.F90 +++ b/components/homme/src/share/physics_mod.F90 @@ -9,7 +9,7 @@ module physics_mod ! ======================= use kinds, only : real_kind ! ======================= - use physical_constants, only : rgas, Rwater_vapor, kappa, g, Rd_on_Rv, Cp, Cpd_on_Cpv, cpwater_vapor + use physical_constants, only : rgas, Rwater_vapor, kappa, gravit, Rd_on_Rv, Cp, Cpd_on_Cpv, cpwater_vapor ! ======================= use physical_constants, only : rearth,p0 ! ======================= diff --git a/components/homme/src/share/prim_advection_base.F90 b/components/homme/src/share/prim_advection_base.F90 index ed31c089d521..a8038ea10623 100644 --- a/components/homme/src/share/prim_advection_base.F90 +++ b/components/homme/src/share/prim_advection_base.F90 @@ -47,7 +47,7 @@ module prim_advection_base ! use kinds, only : real_kind use dimensions_mod, only : nlev, nlevp, np, qsize - use physical_constants, only : rgas, Rwater_vapor, kappa, g, rearth, rrearth, cp + use physical_constants, only : rgas, Rwater_vapor, kappa, gravit, rearth, rrearth, cp use derivative_mod, only : derivative_t, gradient_sphere, divergence_sphere use element_mod, only : element_t use hybvcoord_mod, only : hvcoord_t @@ -229,6 +229,7 @@ subroutine precompute_divdp( elem , hybrid , deriv , dt , nets , nete , n0_qdp ) real(kind=real_kind) , intent(in ) :: dt integer , intent(in ) :: nets , nete , n0_qdp integer :: ie , k + do ie = nets , nete do k = 1 , nlev ! div( U dp Q), @@ -279,6 +280,7 @@ subroutine euler_step( np1_qdp , n0_qdp , dt , elem , hvcoord , hybrid , deriv , limiter_optim_iter_full, limiter_clip_and_sum use bndry_mod , only : bndry_exchangev use hybvcoord_mod , only : hvcoord_t + use deep_atm_mod, only: r_hat_from_phi implicit none integer , intent(in ) :: np1_qdp, n0_qdp real (kind=real_kind), intent(in ) :: dt @@ -299,6 +301,7 @@ subroutine euler_step( np1_qdp , n0_qdp , dt , elem , hvcoord , hybrid , deriv , real(kind=real_kind), dimension(np,np ,nlev ) :: Qtens real(kind=real_kind), dimension(np,np ,nlev ) :: dp,dp_star real(kind=real_kind), dimension(np,np ,nlev,qsize,nets:nete) :: Qtens_biharmonic + real(kind=real_kind), dimension(np,np) :: recip_r_hat real(kind=real_kind), pointer, dimension(:,:,:) :: DSSvar integer :: ie,q,i,j,k, kptr integer :: rhs_viss @@ -476,8 +479,9 @@ subroutine euler_step( np1_qdp , n0_qdp , dt , elem , hvcoord , hybrid , deriv , do q = 1 , qsize do k = 1 , nlev ! dp_star used as temporary instead of divdp (AAM) ! div( U dp Q), - gradQ(:,:,1) = Vstar(:,:,1,k) * elem(ie)%state%Qdp(:,:,k,q,n0_qdp) - gradQ(:,:,2) = Vstar(:,:,2,k) * elem(ie)%state%Qdp(:,:,k,q,n0_qdp) + recip_r_hat = 1.0_real_kind/r_hat_from_phi((elem(ie)%state%phinh_i(:,:,k,n0_qdp) + elem(ie)%state%phinh_i(:,:,k+1,n0_qdp))/2) !DA_CHANGE + gradQ(:,:,1) = Vstar(:,:,1,k) * elem(ie)%state%Qdp(:,:,k,q,n0_qdp) * recip_r_hat + gradQ(:,:,2) = Vstar(:,:,2,k) * elem(ie)%state%Qdp(:,:,k,q,n0_qdp) * recip_r_hat dp_star(:,:,k) = divergence_sphere( gradQ , deriv , elem(ie) ) Qtens(:,:,k) = elem(ie)%state%Qdp(:,:,k,q,n0_qdp) - dt * dp_star(:,:,k) ! optionally add in hyperviscosity computed above: diff --git a/components/homme/src/share/prim_driver_base.F90 b/components/homme/src/share/prim_driver_base.F90 index 27df23dba4d4..8e24dd7cb49c 100644 --- a/components/homme/src/share/prim_driver_base.F90 +++ b/components/homme/src/share/prim_driver_base.F90 @@ -748,6 +748,7 @@ subroutine prim_init2(elem, hybrid, nets, nete, tl, hvcoord) use prim_advection_mod, only: prim_advec_init2 use model_init_mod, only: model_init2 use time_mod, only: timelevel_t, tstep, timelevel_init, nendstep, smooth, nsplit, TimeLevel_Qdp + use deep_atm_mod, only: r_hat_from_phi use control_mod, only: smooth_phis_numcycle #ifdef TRILINOS @@ -775,6 +776,7 @@ subroutine prim_init2(elem, hybrid, nets, nete, tl, hvcoord) real (kind=real_kind) :: dp real (kind=real_kind) :: ps(np,np) ! surface pressure + real (kind=real_kind) :: r_hat(np,np) character(len=80) :: fname character(len=8) :: njusn @@ -929,7 +931,7 @@ end subroutine noxinit ! enddo endif !runtype - + #endif !$OMP MASTER if (runtype==2) then @@ -947,14 +949,12 @@ end subroutine noxinit ! have access to hvcoord to compute dp3d: do ie=nets,nete do k=1,nlev - elem(ie)%state%dp3d(:,:,k,tl%n0)=& - ( hvcoord%hyai(k+1) - hvcoord%hyai(k) )*hvcoord%ps0 + & + elem(ie)%state%dp3d(:,:,k,tl%n0)= ( hvcoord%hyai(k+1) - hvcoord%hyai(k) )*hvcoord%ps0 + & ( hvcoord%hybi(k+1) - hvcoord%hybi(k) )*elem(ie)%state%ps_v(:,:,tl%n0) enddo end do #endif - ! For new runs, and branch runs, convert state variable Q to (Qdp) ! because initial conditon reads in Q, not Qdp ! restart runs will read dpQ from restart file @@ -978,7 +978,6 @@ end subroutine noxinit enddo enddo endif - call model_init2(elem(:), hybrid,deriv1,hvcoord,tl,nets,nete) ! advective and viscious CFL estimates @@ -1060,6 +1059,7 @@ subroutine prim_run_subcycle(elem, hybrid,nets,nete, dt, single_column, tl, hvco use hybvcoord_mod, only: hvcoord_t use parallel_mod, only: abortmp use prim_state_mod, only: prim_printstate +use eos, only: pnh_and_exner_from_eos use vertremap_mod, only: vertical_remap use reduction_mod, only: parallelmax use time_mod, only: TimeLevel_t, timelevel_update, timelevel_qdp, nsplit, tstep @@ -1083,6 +1083,7 @@ subroutine prim_run_subcycle(elem, hybrid,nets,nete, dt, single_column, tl, hvco real(kind=real_kind) :: dp_np1(np,np) integer :: ie,i,j,k,n,q,t,scm_dum integer :: n0_qdp,np1_qdp,r,nstep_end,nets_in,nete_in,step_factor + real (kind=real_kind) :: pnh(np, np, nlev), exner(np, np, nlev), dpnh_i(np,np,nlevp) logical :: compute_diagnostics ! compute timesteps for tracer transport and vertical remap @@ -1107,7 +1108,8 @@ subroutine prim_run_subcycle(elem, hybrid,nets,nete, dt, single_column, tl, hvco ! compute scalar diagnostics if currently active if (compute_diagnostics) call run_diagnostics(elem,hvcoord,tl,3,.true.,nets,nete) - if (prim_step_type == 1) then + + if (prim_step_type == 1) then call TimeLevel_Qdp(tl, dt_tracer_factor, n0_qdp, np1_qdp) #if !defined(CAM) && !defined(SCREAM) @@ -1116,8 +1118,8 @@ subroutine prim_run_subcycle(elem, hybrid,nets,nete, dt, single_column, tl, hvco ! homme. call compute_test_forcing(elem,hybrid,hvcoord,tl%n0,n0_qdp,dt_remap,nets,nete,tl) #endif - - call applyCAMforcing_remap(elem,hvcoord,tl%n0,n0_qdp,dt_remap,nets,nete) + + !call applyCAMforcing_remap(elem,hvcoord,tl%n0,n0_qdp,dt_remap,nets,nete) ! E(1) Energy after CAM forcing if (compute_diagnostics) call run_diagnostics(elem,hvcoord,tl,1,.true.,nets,nete) @@ -1175,7 +1177,7 @@ subroutine prim_run_subcycle(elem, hybrid,nets,nete, dt, single_column, tl, hvco nets_in=nets nete_in=nete endif - + call vertical_remap(hybrid,elem,hvcoord,dt_remap,tl%np1,np1_qdp,nets_in,nete_in) elseif(prim_step_type == 2) then ! This time stepping routine permits the vertical remap time @@ -1257,7 +1259,6 @@ subroutine prim_step(elem, hybrid,nets,nete, dt, tl, hvcoord, compute_diagnostic real (kind=real_kind) :: maxcflx, maxcfly real (kind=real_kind) :: dp_np1(np,np) logical :: compute_diagnostics - dt_q = dt*dt_tracer_factor call set_tracer_transport_derived_values(elem, nets, nete, tl) @@ -1322,6 +1323,8 @@ subroutine prim_step_flexible(hybrid, elem, nets, nete, dt, tl, hvcoord, compute use prim_state_mod, only: prim_printstate use vertremap_mod, only: vertical_remap use sl_advection, only: sl_vertically_remap_tracers + use physical_constants, only: rearth, gravit + use deep_atm_mod, only: r_hat_from_phi type(element_t), intent(inout) :: elem(:) type(hybrid_t), intent(in) :: hybrid ! distributed parallel structure (shared) @@ -1333,6 +1336,7 @@ subroutine prim_step_flexible(hybrid, elem, nets, nete, dt, tl, hvcoord, compute logical, intent(in) :: compute_diagnostics real(kind=real_kind) :: dt_q, dt_remap, dp(np,np,nlev) + real(kind=real_kind) :: r_hat(np,np) integer :: ie, q, k, n, n0_qdp, np1_qdp logical :: compute_diagnostics_it, apply_forcing @@ -1404,7 +1408,7 @@ subroutine prim_step_flexible(hybrid, elem, nets, nete, dt, tl, hvcoord, compute ! Prescribed winds are evaluated on reference levels, ! not floating levels, so don't remap, just update dp3d. do ie = nets,nete - elem(ie)%state%ps_v(:,:,tl%np1) = hvcoord%hyai(1)*hvcoord%ps0 + & + elem(ie)%state%ps_v(:,:,tl%np1) = hvcoord%hyai(1)*hvcoord%ps0 + & sum(elem(ie)%state%dp3d(:,:,:,tl%np1),3) do k=1,nlev dp(:,:,k) = (hvcoord%hyai(k+1) - hvcoord%hyai(k))*hvcoord%ps0 + & @@ -1569,9 +1573,10 @@ subroutine applyCAMforcing_tracers(elem,hvcoord,np1,np1_qdp,dt,adjustment) ! use control_mod, only : use_moisture, dt_remap_factor use hybvcoord_mod, only : hvcoord_t + use deep_atm_mod, only: r_hat_from_phi, quasi_hydrostatic_terms #ifdef MODEL_THETA_L use control_mod, only : theta_hydrostatic_mode - use physical_constants, only : cp, g, kappa, Rgas, p0 + use physical_constants, only : cp, gravit, kappa, Rgas, p0, rearth use element_ops, only : get_temperature, get_r_star, get_hydro_pressure use eos, only : pnh_and_exner_from_eos #ifdef HOMMEXX_BFB_TESTING @@ -1600,6 +1605,7 @@ subroutine applyCAMforcing_tracers(elem,hvcoord,np1,np1_qdp,dt,adjustment) real (kind=real_kind) :: rstarn1(np,np,nlev) real (kind=real_kind) :: exner(np,np,nlev) real (kind=real_kind) :: dpnh_dp_i(np,np,nlevp) + real (kind=real_kind) :: r_hat(np,np) #endif #ifdef HOMMEXX_BFB_TESTING @@ -1720,8 +1726,8 @@ subroutine applyCAMforcing_tracers(elem,hvcoord,np1,np1_qdp,dt,adjustment) if (adjust_ps) then ! compute new dp3d from adjusted ps() do k=1,nlev - dp_adj(:,:,k) = ( hvcoord%hyai(k+1) - hvcoord%hyai(k) )*hvcoord%ps0 + & - ( hvcoord%hybi(k+1) - hvcoord%hybi(k))*ps(:,:) + dp_adj(:,:,k) = (( hvcoord%hyai(k+1) - hvcoord%hyai(k) )*hvcoord%ps0 + & + ( hvcoord%hybi(k+1) - hvcoord%hybi(k))*ps(:,:)) enddo endif elem%state%dp3d(:,:,:,np1)=dp_adj(:,:,:) diff --git a/components/homme/src/share/prim_si_mod.F90 b/components/homme/src/share/prim_si_mod.F90 index 9e2c8bec9a26..f309504804af 100644 --- a/components/homme/src/share/prim_si_mod.F90 +++ b/components/homme/src/share/prim_si_mod.F90 @@ -567,7 +567,7 @@ end subroutine geopotential_t subroutine prim_set_mass(elem, tl,hybrid,hvcoord,nets,nete) use kinds, only : real_kind use control_mod, only : initial_total_mass - use physical_constants, only : g + use physical_constants, only : gravit use element_mod, only : element_t use time_mod, only : timelevel_t use hybvcoord_mod, only : hvcoord_t @@ -592,7 +592,7 @@ subroutine prim_set_mass(elem, tl,hybrid,hvcoord,nets,nete) nm1=tl%nm1 np1=tl%np1 - scale=1/g ! assume code is using Pa + scale=1/gravit ! assume code is using Pa Warning: not adapted to Deep Atmosphere if (hvcoord%ps0 < 2000 ) scale=100*scale ! code is using mb ! after scaling, Energy is in J/m**2, Mass kg/m**2 diff --git a/components/homme/src/test_src/.dcmip2016-supercell.F90.swp b/components/homme/src/test_src/.dcmip2016-supercell.F90.swp new file mode 100644 index 000000000000..f048c2486190 Binary files /dev/null and b/components/homme/src/test_src/.dcmip2016-supercell.F90.swp differ diff --git a/components/homme/src/test_src/asp_tests.F90 b/components/homme/src/test_src/asp_tests.F90 index aac88e578fbc..18e963ffb8f8 100644 --- a/components/homme/src/test_src/asp_tests.F90 +++ b/components/homme/src/test_src/asp_tests.F90 @@ -982,7 +982,7 @@ module asp_tests use jw, only: u_wind, v_wind, temperature, surface_geopotential, tracer_q1_q2,& tracer_q3, perturbation_longitude, perturbation_latitude, deg2rad ! _EXTERNAL use parallel_mod, only: abortmp -use physical_constants, only: p0, g +use physical_constants, only: p0, g=>gravit implicit none diff --git a/components/homme/src/test_src/baroclinic_inst_mod.F90 b/components/homme/src/test_src/baroclinic_inst_mod.F90 index 32c3381ef524..67d2cf6efeb9 100644 --- a/components/homme/src/test_src/baroclinic_inst_mod.F90 +++ b/components/homme/src/test_src/baroclinic_inst_mod.F90 @@ -22,7 +22,7 @@ module baroclinic_inst_mod ! ==================== use kinds, only : real_kind, iulog ! ==================== - use physical_constants, only : omega, rearth, rgas, p0, dd_pi, Cp,g + use physical_constants, only : omega, rearth, rgas, p0, dd_pi, Cp,g=>gravit ! ==================== use dimensions_mod, only : nlev,np, qsize ! ==================== diff --git a/components/homme/src/test_src/dcmip12_wrapper.F90 b/components/homme/src/test_src/dcmip12_wrapper.F90 index bd4da6eff78a..c8152b62b638 100644 --- a/components/homme/src/test_src/dcmip12_wrapper.F90 +++ b/components/homme/src/test_src/dcmip12_wrapper.F90 @@ -91,7 +91,7 @@ subroutine dcmip2012_test1_1(elem,hybrid,hvcoord,nets,nete,time,n0,n1) dp = pressure_thickness(ps,k,hvcoord) call set_state(u,v,w,T,ps,phis,p,dp,zm(k),g, i,j,k,elem(ie),n0,n1) - if(time==0) call set_tracers(q,qsize,dp,i,j,k,lat,lon,elem(ie)) + if(time==0) call set_tracers(q,qsize,dp,zm(k),i,j,k,lat,lon,elem(ie)) enddo; enddo; enddo; enddo @@ -167,7 +167,7 @@ subroutine dcmip2012_test1_1_conv(elem,hybrid,hvcoord,nets,nete,time,n0,n1) dp = pressure_thickness(ps,k,hvcoord) call set_state(u,v,w,T,ps,phis,p,dp,zm(k),g, i,j,k,elem(ie),n0,n1) - if(time==0) call set_tracers(q,qsize,dp,i,j,k,lat,lon,elem(ie)) + if(time==0) call set_tracers(q,qsize,dp,zm(k),i,j,k,lat,lon,elem(ie)) enddo; enddo; enddo; enddo @@ -235,7 +235,7 @@ subroutine dcmip2012_test1_2(elem,hybrid,hvcoord,nets,nete,time,n0,n1) call test1_advection_hadley(time,lon,lat,p,z,zcoords,u,v,w,t,phis,ps,rho,q(1),q(2)) dp = pressure_thickness(ps,k,hvcoord) call set_state(u,v,w,T,ps,phis,p,dp,zm(k),g, i,j,k,elem(ie),n0,n1) - if(time==0) call set_tracers(q,qsize,dp,i,j,k,lat,lon,elem(ie)) + if(time==0) call set_tracers(q,qsize,dp,zm(k),i,j,k,lat,lon,elem(ie)) enddo; enddo; enddo; enddo @@ -306,7 +306,7 @@ subroutine dcmip2012_test1_3(elem,hybrid,hvcoord,nets,nete,time,n0,n1,deriv) call test1_advection_orography(lon,lat,p,z,zcoords,cfv,use_eta,hyam,hybm,gc,u,v,w,t,phis,ps,rho,q(1),q(2),q(3),q(4)) dp = pressure_thickness(ps,k,hvcoord) call set_state(u,v,w,T,ps,phis,p,dp,zm(k),g, i,j,k,elem(ie),n0,n1) - if(time==0) call set_tracers(q,qsize,dp,i,j,k,lat,lon,elem(ie)) + if(time==0) call set_tracers(q,qsize,dp,zm(k),i,j,k,lat,lon,elem(ie)) enddo; enddo; enddo; enddo @@ -373,7 +373,7 @@ subroutine dcmip2012_test2_0(elem,hybrid,hvcoord,nets,nete) !let's get an analytical \phi he = (T0 - T)/gamma call set_state(u,v,w,T,ps,phis,p,dp,he,g, i,j,k,elem(ie),1,nt) - call set_tracers(q,qsize,dp,i,j,k,lat,lon,elem(ie)) + call set_tracers(q,qsize,dp,he,i,j,k,lat,lon,elem(ie)) enddo; enddo; enddo; do k=1,nlevp; do j=1,np; do i=1,np call get_coordinates(lat,lon,hyai,hybi, i,j,k,elem(ie),hvcoord) @@ -432,7 +432,7 @@ subroutine dcmip2012_test2_x(elem,hybrid,hvcoord,nets,nete,shear) ! call set_state(u,v,w,T,ps,phis,p,dp,zm(k),g, i,j,k,elem(ie),1,nt) ! This test obtains analytical height and returns it, so, we use it for \phi ... call set_state(u,v,w,T,ps,phis,p,dp,z,g, i,j,k,elem(ie),1,nt) - call set_tracers(q,qsize,dp,i,j,k,lat,lon,elem(ie)) + call set_tracers(q,qsize,dp,z,i,j,k,lat,lon,elem(ie)) ! ... or we can use discrete hydro state to init \phi. enddo; enddo; enddo; @@ -602,7 +602,7 @@ subroutine dcmip2012_test3(elem,hybrid,hvcoord,nets,nete) call test3_gravity_wave(lon,lat,p,z,zcoords,use_eta,hyam,hybm,u,v,w,T,T_mean,phis,ps,rho,rho_mean,q(1)) dp = pressure_thickness(ps,k,hvcoord) call set_state(u,v,w,T,ps,phis,p,dp,zm(k),g, i,j,k,elem(ie),1,nt) - call set_tracers(q,qsize, dp,i,j,k,lat,lon,elem(ie)) + call set_tracers(q,qsize, dp,z,i,j,k,lat,lon,elem(ie)) enddo; enddo; enddo; do k=1,nlevp; do j=1,np; do i=1,np call get_coordinates(lat,lon,hyai,hybi, i,j,k,elem(ie),hvcoord) @@ -662,7 +662,7 @@ subroutine dcmip2012_test4_init(elem,hybrid,hvcoord,nets,nete) !init only <=qsize tracers qs = min(qsize,3) - call set_tracers(qarray(1:qs),qs, dp,i,j,k,lat,lon,elem(ie)) + call set_tracers(qarray(1:qs),qs, dp,z,i,j,k,lat,lon,elem(ie)) enddo; enddo; enddo; do k=1,nlevp @@ -785,30 +785,34 @@ real(rl) function pressure_thickness(ps,k,hv) !_____________________________________________________________________ -subroutine set_tracers(q,nq, dp,i,j,k,lat,lon,elem) - +subroutine set_tracers(q,nq, dp,zm,i,j,k,lat,lon,elem) + use physical_constants, only: rearth + use deep_atm_mod, only: r_hat_from_phi, phi_from_z ! set tracer values at node(i,j,k) - real(rl), intent(in) :: q(nq), dp, lat, lon + real(rl), intent(in) :: q(nq), dp,zm, lat, lon integer, intent(in) :: i,j,k,nq type(element_t), intent(inout) :: elem - + real(rl) :: r_hat, pdensity real(rl), parameter :: wl = 1.0 ! checkerboard wavelength in dg integer :: qi + r_hat = r_hat_from_phi(phi_from_z(zm)) + pdensity = dp + if (nq>qsize) call abortmp('qsize set too small for dcmip test case') ! set known tracers to q and the rest to a checkerboard pattern elem%state%Q(i,j,k,1:nq) = q ! compute tracer mass qdp from mixing ratio q do qi = 1,nq - elem%state%Qdp (i,j,k,qi,:) = q(qi)*dp + elem%state%Qdp (i,j,k,qi,:) = q(qi)*pdensity enddo ! set any remaining tracers to 1 do qi = nq+1,qsize elem%state%Q(i,j,k,qi) = 1 - elem%state%Qdp (i,j,k,qi,:) = elem%state%Q(i,j,k,qi)*dp + elem%state%Qdp (i,j,k,qi,:) = elem%state%Q(i,j,k,qi)*pdensity enddo diff --git a/components/homme/src/test_src/dcmip16_wrapper.F90 b/components/homme/src/test_src/dcmip16_wrapper.F90 index 719710e0bb27..ddefd9a8f7ff 100644 --- a/components/homme/src/test_src/dcmip16_wrapper.F90 +++ b/components/homme/src/test_src/dcmip16_wrapper.F90 @@ -25,7 +25,8 @@ module dcmip16_wrapper use parallel_mod, only: abortmp,iam use element_ops, only: set_state, set_state_i, set_elem_state, get_state, tests_finalize,& set_forcing_rayleigh_friction -use physical_constants, only: p0, g, Rgas, kappa, Cp, Rwater_vapor, pi=>dd_pi +use physical_constants, only: p0, g=>gravit, Rgas, kappa, Cp, Rwater_vapor, pi=>dd_pi, rearth, omega +use deep_atm_mod, only: atm_is_deep use reduction_mod, only: parallelmax, parallelmin use terminator, only: initial_value_terminator, tendency_terminator use time_mod, only: time_at, TimeLevel_t @@ -42,7 +43,8 @@ module dcmip16_wrapper real(rl) :: sample_period = 2.0_rl real(rl) :: rad2dg = 180.0_rl/pi - +real(rl), parameter :: gam = 0.005 +real(rl), parameter :: T0 = 300 type :: PhysgridData_t integer :: nphys real(rl), allocatable :: ps(:,:), zs(:,:), T(:,:,:), uv(:,:,:,:), omega_p(:,:,:), q(:,:,:,:) @@ -71,7 +73,8 @@ end subroutine dcmip2016_init !_____________________________________________________________________ subroutine dcmip2016_test1(elem,hybrid,hvcoord,nets,nete) - + use deep_atm_mod, only: r_hat_from_phi, phi_from_z, z_from_phi, g_from_phi + use eos, only: pnh_and_exner_from_eos ! moist baroclinic wave test type(element_t), intent(inout), target :: elem(:) ! element array @@ -80,19 +83,28 @@ subroutine dcmip2016_test1(elem,hybrid,hvcoord,nets,nete) integer, intent(in) :: nets,nete ! start, end element index integer, parameter :: use_zcoords = 0 ! use vertical pressure coordinates - integer, parameter :: is_deep = 0 ! use shallow atmosphere approximation + logical, parameter :: is_deep = atm_is_deep ! use shallow atmosphere approximation integer, parameter :: pertt = 0 ! use exponential perturbation type real(rl), parameter :: dcmip_X = 1.0_rl ! full scale planet integer :: moist ! use moist version - integer :: i,j,k,ie ! loop indices + integer :: i,j,k,ie, newton_idx ! loop indices real(rl):: lon,lat ! pointwise coordiantes - real(rl), dimension(np,np,nlev):: p,z,u,v,w,T,thetav,rho,dp ! field values - real(rl), dimension(np,np,nlevp):: p_i,w_i,z_i + real(rl), dimension(np,np,nlev):: p,dptmpp,z,u,v,w,T,thetav,rho,dp ! field values + real(rl), dimension(np,np,nlev):: dphi_tmp ! field values + real(rl), dimension(np,np,nlevp):: p_i,ptmp_i,w_i,z_i,mu + real(rl), dimension(np, np,nlev) :: phi_guess, residual, rhs, deriv, dmass, dmass_i, r_hat, r_hat_pert, & + pnh_next, pnh_next_pert, dphi_hack, pnh_lower, z_guess, z_guess_pert, t_pert, pnh,exner, u_pert + real(rl), dimension(np,np,nlevp) :: phi_prev, pnh_prev, dpnh_dp_i real(rl), dimension(np,np):: ps, phis + real(rl) :: u_tmp, v_tmp, q_tmp, p_tmp, t_tmp real(rl), dimension(np,np,nlev,6):: q + integer, parameter :: max_newton = 400 + real(rl), parameter :: newton_eps = 1e-8 real(rl) :: min_thetav, max_thetav + logical :: simple_atm_hack = .false., do_jank_rootfinding=.true. + logical, parameter :: mu_instead_of_eos_rootfinding = .false. min_thetav = +huge(rl) max_thetav = -huge(rl) @@ -102,55 +114,197 @@ subroutine dcmip2016_test1(elem,hybrid,hvcoord,nets,nete) if (hybrid%masterthread) write(iulog,*) 'initializing dcmip2016 test 1: moist baroclinic wave' if (qsize<6) call abortmp('ERROR: test requires qsize>=6') - - ! set initial conditions + do ie = nets,nete - do k=1,nlevp; do j=1,np; do i=1,np - ! no surface topography - p_i(i,j,k) = p0*hvcoord%etai(k) + do k=1,nlev + dp(:,:, k) =-( (hvcoord%etai(k)-hvcoord%etai(k+1)) * p0) + end do + do k=1,nlevp; do j=1,np; do i=1,np lon = elem(ie)%spherep(i,j)%lon lat = elem(ie)%spherep(i,j)%lat - w_i(i,j,k) = 0.0d0 - ! call this only to compute z_i, will ignore all other output - call baroclinic_wave_test(is_deep,moist,pertt,dcmip_X,lon,lat,p_i(i,j,k),& + + p_i(i,j,k) = p0*hvcoord%etai(k) + call baroclinic_wave_test(is_deep,moist,pertt,dcmip_X,lon,lat,p_i(i,j,k),& z_i(i,j,k),use_zcoords,u(i,j,1),v(i,j,1),T(i,j,1),thetav(i,j,1),phis(i,j),ps(i,j),rho(i,j,1),q(i,j,1,1)) - enddo; enddo; enddo - do k=1,nlev; do j=1,np; do i=1,np + if (k > 1) then + dphi_tmp(i,j,k-1) = (phi_from_z(z_i(i,j,k-1)) - phi_from_z(z_i(i,j,k)))/r_hat_from_phi(phi_from_z(z_i(i,j,k)))**2 + end if + end do; end do; end do + q(:,:,:,1:5) = 0.0d0 + q(:,:,:,6) = 1 + w(:,:,:) = 0.0d0 + w_i(:,:,:) = 0.0d0 + + + dphi_tmp = 10 + + phi_prev(:,:,nlevp) = phis + z_i(:, :, nlevp) = z_from_phi(phis) + pnh_prev(:,:,nlevp) = ps + do k=nlev,1,-1; do j=1,np; do i=1,np + + lon = elem(ie)%spherep(i,j)%lon + lat = elem(ie)%spherep(i,j)%lat + + +! =================================== + if (k.eq. nlev ) then + dmass_i(i,j,k) = -dp(i,j,k)/2 + else + dmass_i(i,j,k) = -(dp(i,j,k)+dp(i,j,k+1))/2 + end if + dmass(i,j,k) = -dp(i,j,k) + + phi_guess(i,j,k) = phi_prev(i,j,k+1) + dphi_tmp(i,j,k) + + if (mu_instead_of_eos_rootfinding) then + + do newton_idx=1, max_newton + z_guess(i,j,k) = z_from_phi((phi_prev(i,j,k+1) + phi_guess(i,j,k))/2) + z_guess_pert(i,j,k) = z_from_phi((phi_prev(i,j,k+1) +phi_guess(i,j,k)+newton_eps)/2) + call baroclinic_wave_test(is_deep,moist,pertt,dcmip_X,lon,lat, & + pnh_next(i,j,k),z_guess(i,j,k),1,u(i,j,k),v(i,j,k),T(i,j,k),thetav(i,j,k),phis(i,j),ps(i,j),rho(i,j,k),q(i,j,k,1)) + call baroclinic_wave_test(is_deep,moist,pertt,dcmip_X,lon,lat, & + pnh_next_pert(i,j,k),z_guess_pert(i,j,k),1,u_pert(i,j,k),v(i,j,k),t_pert(i,j,k),thetav(i,j,k),phis(i,j),ps(i,j),rho(i,j,k),q(i,j,k,1)) + + + r_hat(i,j,k) = r_hat_from_phi(( phi_prev(i,j,k+1) + phi_guess(i,j,k))/2) + r_hat_pert(i,j,k) = r_hat_from_phi(( phi_prev(i,j,k+1) + phi_guess(i,j,k)+newton_eps)/2) + + pnh_next(i,j,k) = -Rgas * T(i,j,k) * dmass(i,j,k) / (r_hat(i,j,k)**2 * (phi_guess(i,j,k) - phi_prev(i,j,k+1))) + pnh_next_pert(i,j,k) = -Rgas * t_pert(i,j,k) * dmass(i,j,k) / (r_hat_pert(i,j,k)**2 * (phi_guess(i,j,k) + newton_eps - phi_prev(i,j,k+1))) + !end if + residual(i,j,k) = -u(i,j,k)**2/(z_guess(i,j,k)+rearth) - 2 * omega * cos(lat) * u(i,j,k) + & + g_from_phi(phi_guess(i,j,k)) * (1 - (pnh_next(i,j,k)-pnh_prev(i,j,k+1))/ dmass_i(i,j,k) * r_hat_from_phi(phi_prev(i,j,k+1))**2 ) + if (abs(residual(i,j,k)) < 1e-10) then + exit + end if + deriv(i,j,k) = ( (-u_pert(i,j,k)**2/(z_guess_pert(i,j,k)+rearth) - 2 * omega * cos(lat) * u_pert(i,j,k) + & + g_from_phi(phi_guess(i,j,k)) * (1 - (pnh_next_pert(i,j,k)-pnh_prev(i,j,k+1))/ dmass_i(i,j,k) * r_hat_from_phi(phi_prev(i,j,k+1))**2 ))-& + (-u(i,j,k)**2/(z_guess(i,j,k)+rearth) - 2 * omega * cos(lat) * u(i,j,k) + & + g_from_phi(phi_guess(i,j,k)) * (1 - (pnh_next(i,j,k)-pnh_prev(i,j,k+1))/ dmass_i(i,j,k) * r_hat_from_phi(phi_prev(i,j,k+1))**2 )))/newton_eps + + + + if (isnan( residual(i,j,k)/deriv(i,j,k)) .or. phi_guess(i,j,k) < phi_prev(i,j,k+1)) then ! .or. pnh_next(i,j,k) > pnh_prev(i,j,k+1)) then + write(iulog,*)'WARNING: Deep atm geopotential rootfinding encountered zero at lev: ', k + write(iulog, *) "phi next: ", phi_guess(i,j,k), "phi_prev: ", phi_prev(i,j,k+1) + write(iulog, *) "pnh_next: ", pnh_next(i,j,k), "pnh_prev: ", pnh_prev(i,j,k+1) + call abort + end if + + phi_guess(i,j,k) = phi_guess(i,j,k) - residual(i,j,k)/deriv(i,j,k) + - ! no surface topography - p(i,j,k) = p0*hvcoord%etam(k) - dp(i,j,k) = (hvcoord%etai(k+1)-hvcoord%etai(k))*p0 + end do + + else + do newton_idx=1, max_newton + z_guess(i,j,k) = z_from_phi((phi_guess(i,j,k) + phi_prev(i,j,k+1))/2) + z_guess_pert(i,j,k) = z_from_phi((phi_prev(i,j,k+1)+phi_guess(i,j,k)+newton_eps)/2) + call baroclinic_wave_test(is_deep,moist,pertt,dcmip_X,lon,lat, & + p_tmp,z_guess(i,j,k),1,u(i,j,k),v(i,j,k),T(i,j,k),thetav(i,j,k),phis(i,j),ps(i,j),rho(i,j,k),q(i,j,k,1)) + call baroclinic_wave_test(is_deep,moist,pertt,dcmip_X,lon,lat, & + pnh_next_pert(i,j,k),z_guess_pert(i,j,k),1,u(i,j,k),v(i,j,k),t_pert(i,j,k),thetav(i,j,k),phis(i,j),ps(i,j),rho(i,j,k),q(i,j,k,1)) + + r_hat(i,j,k) = r_hat_from_phi((phi_guess(i,j,k) + phi_prev(i,j,k+1))/2) + r_hat_pert(i,j,k) = r_hat_from_phi((phi_prev(i,j,k+1) + phi_guess(i,j,k)+newton_eps)/2) - lon = elem(ie)%spherep(i,j)%lon - lat = elem(ie)%spherep(i,j)%lat - q(i,j,k,1:5) = 0.0d0 - q(i,j,k,6) = 1 - w(i,j,k) = 0.0d0 - call baroclinic_wave_test(is_deep,moist,pertt,dcmip_X,lon,lat,p(i,j,k),& - z(i,j,k),use_zcoords,u(i,j,k),v(i,j,k),T(i,j,k),thetav(i,j,k),phis(i,j),ps(i,j),rho(i,j,k),q(i,j,k,1)) + pnh_next(i,j,k) = -Rgas * T(i,j,k) * dmass(i,j,k) / (r_hat(i,j,k)**2 * (phi_guess(i,j,k) - phi_prev(i,j,k+1))) + pnh_next_pert(i,j,k) = -Rgas * t_pert(i,j,k) * dmass(i,j,k) / (r_hat_pert(i,j,k)**2 * (phi_guess(i,j,k) + newton_eps - phi_prev(i,j,k+1))) + residual(i,j,k) = pnh_next(i,j,k) - p_tmp + if (abs(residual(i,j,k)) < 1e-10) then + exit + end if + deriv(i,j,k) = (pnh_next_pert(i,j,k) & + - pnh_next(i,j,k))/newton_eps - ! initialize tracer chemistry - call initial_value_terminator( lat*rad2dg, lon*rad2dg, q(i,j,k,4), q(i,j,k,5) ) - call set_tracers(q(i,j,k,1:6),6,dp(i,j,k),i,j,k,lat,lon,elem(ie)) + if (isnan( residual(i,j,k)/deriv(i,j,k)) .or. phi_guess(i,j,k) < phi_prev(i,j,k+1)) then + write(iulog,*)'WARNING: Deep atm geopotential rootfinding encountered zero at lev: ', k + write(iulog, *) "phi next: ", phi_guess(i,j,k), "phi_prev: ", phi_prev(i,j,k+1) + write(iulog, *) "pnh_next: ", pnh_next(i,j,k), "pnh_prev: ", pnh_prev(i,j,k+1) + call abort + end if - min_thetav = min( min_thetav, thetav(i,j,k) ) - max_thetav = max( max_thetav, thetav(i,j,k) ) + phi_guess(i,j,k) = phi_guess(i,j,k) - residual(i,j,k)/deriv(i,j,k) - enddo; enddo; enddo + end do + + + + end if + if (abs(residual(i,j,k)) > 1e-9) then + write(*,*)'WARNING: Deep atm geopotential rootfinding failed at level ', k, ' residual: ', abs(residual(i,j,k)) + call abort + end if + pnh_prev(i,j,k) = pnh_next(i,j,k) + p(i,j,k) = pnh_next(i,j,k) + z_i(i, j, k) = z_from_phi(phi_guess(i, j,k)) + z(i,j,k) = z_from_phi((phi_guess(i,j,k) + phi_prev(i,j,k+1))/2) + + phi_prev(i,j,k) = phi_guess(i,j,k) +! =================================== + + end do; end do; end do + !do k=1,nlev + ! write(*,*) "k: ", k, " max pnh: ", maxval(p(:,:,k)), "min pnh: ", minval(p(:,:,k)) + !end do + + +! do k=1,nlev; do j=1,np; do i=1,np +! ! initialize tracer chemistry +! call initial_value_terminator( lat*rad2dg, lon*rad2dg, q(i,j,k,4), q(i,j,k,5) ) +! call set_tracers(q(i,j,k,1:6),6,dp(i,j,k),z(i,j,k),i,j,k,lat,lon,elem(ie)) +! +! min_thetav = min( min_thetav, thetav(i,j,k) ) +! max_thetav = max( max_thetav, thetav(i,j,k) ) +! +! enddo; enddo; enddo + call set_elem_state(u,v,w,w_i,T,ps,phis,p,dp,z,z_i,g,elem(ie),1,nt,ntQ=1) call tests_finalize(elem(ie),hvcoord) + + call pnh_and_exner_from_eos(hvcoord,elem(ie)%state%vtheta_dp(:,:,:,nt),& + elem(ie)%state%dp3d(:,:,:,nt),elem(ie)%state%phinh_i(:,:,:,nt),pnh,exner,dpnh_dp_i) +! call pnh_and_exner_from_eos(hvcoord,T*dp*((p/p0)**(-kappa)),& +! dp,phi_from_z(z_i, nlevp),pnh,exner,dpnh_dp_i) + + do k=1,nlevp + if (maxval(abs(dpnh_dp_i(:,:,k)-1)) > 1e-9) then + write(*,*) "dpnh_dp_i error at lev ", k, ": ", maxval(dpnh_dp_i(:,:,k)-1), "minval: ", minval(dpnh_dp_i(:,:,k)-1) + !call abort + end if + end do enddo sample_period = 1800.0 ! sec !print *,"min thetav = ",min_thetav, "max thetav=",max_thetav +end subroutine +subroutine z_given_p(p, z) +real(rl), intent(in) :: p +real(rl), intent(out) :: z + if (gam > 1e-6) then + z = T0/max(1e-6, gam) * (1.0 - (p/p0)**(max(1e-6, gam)*Rgas/g)) + else + z = Rgas/g * T0* log(p0/p) + end if +end subroutine +subroutine t_given_z(z, t) +real(rl), intent(in) :: z +real(rl), intent(out) :: t + if (gam > 1e-6) then + t= T0 - gam * z + else + t = T0 + end if end subroutine subroutine dcmip2016_test1_pg(elem,hybrid,hvcoord,nets,nete,nphys) @@ -238,7 +392,7 @@ subroutine dcmip2016_test2(elem,hybrid,hvcoord,nets,nete) q(2)=0 q(3)=0 - call set_tracers(q(:),3,dp(i,j,k),i,j,k,lat,lon,elem(ie)) + call set_tracers(q(:),3,dp(i,j,k),z(i,j,k),i,j,k,lat,lon,elem(ie)) enddo; enddo; enddo; call set_elem_state(u,v,w,w_i,T,ps,phis,p,dp,z,z_i,g,elem(ie),1,nt,ntQ=1) @@ -344,7 +498,7 @@ subroutine dcmip2016_test3(elem,hybrid,hvcoord,nets,nete) q (i,j,k,2)= 0 ! no initial clouds q (i,j,k,3)= 0 ! no initial rain - call set_tracers(q(i,j,k,:),3,dp(i,j,k),i,j,k,lat,lon,elem(ie)) + call set_tracers(q(i,j,k,:),3,dp(i,j,k),z(i,j,k),i,j,k,lat,lon,elem(ie)) enddo; enddo enddo diff --git a/components/homme/src/test_src/dcmip2012_test1_conv.F90 b/components/homme/src/test_src/dcmip2012_test1_conv.F90 index 274dfcacb588..e50ea46bb189 100644 --- a/components/homme/src/test_src/dcmip2012_test1_conv.F90 +++ b/components/homme/src/test_src/dcmip2012_test1_conv.F90 @@ -9,7 +9,7 @@ SUBROUTINE test1_conv_advection_deformation (time,lon,lat,p,z,zcoords,u,v,w,t,ph !----------------------------------------------------------------------- ! Use physical constants consistent with HOMME - use physical_constants, only: a=>rearth0, Rd => Rgas, g, cp, pi=>dd_pi, p0 + use physical_constants, only: a=>rearth0, Rd => Rgas, g=>gravit, cp, pi=>dd_pi, p0 real(8), intent(in) :: time ! simulation time (s) real(8), intent(in) :: lon ! Longitude (radians) diff --git a/components/homme/src/test_src/dcmip2016-baroclinic.F90 b/components/homme/src/test_src/dcmip2016-baroclinic.F90 index db02714c0042..db99434bc7f9 100644 --- a/components/homme/src/test_src/dcmip2016-baroclinic.F90 +++ b/components/homme/src/test_src/dcmip2016-baroclinic.F90 @@ -40,7 +40,7 @@ MODULE baroclinic_wave ! !======================================================================= -use physical_constants, only: g0=>g,kappa0=>kappa,Rgas,Cp0=>Cp,Rwater_vapor,rearth0,omega0, dd_pi +use physical_constants, only: g0=>gravit,kappa0=>kappa,Rgas,Cp0=>Cp,Rwater_vapor,rearth0,omega0, dd_pi IMPLICIT NONE @@ -88,9 +88,9 @@ MODULE baroclinic_wave lapse = 0.005d0 ! lapse rate parameter REAL(8), PARAMETER :: & - pertu0 = 0.5d0 , & ! SF Perturbation wind velocity (m/s) + pertu0 = 0.0d0 , & ! SF Perturbation wind velocity (m/s) pertr = 1.d0/6.d0 , & ! SF Perturbation radius (Earth radii) - pertup = 1.0d0 , & ! Exp. perturbation wind velocity (m/s) + pertup = 0.0d0 , & ! Exp. perturbation wind velocity (m/s) pertexpr = 0.1d0 , & ! Exp. perturbation radius (Earth radii) pertlon = pi/9.d0 , & ! Perturbation longitude pertlat = 2.d0*pi/9.d0, & ! Perturbation latitude @@ -122,9 +122,9 @@ SUBROUTINE baroclinic_wave_test(deep,moist,pertt,X,lon,lat,p,z,zcoords,u,v,t,the ! input/output params parameters at given location !----------------------------------------------------------------------- INTEGER, INTENT(IN) :: & - deep, & ! Deep (1) or Shallow (0) test case moist, & ! Moist (1) or Dry (0) test case pertt ! Perturbation type + LOGICAL, INTENT(IN) :: deep REAL(8), INTENT(IN) :: & lon, & ! Longitude (radians) @@ -182,7 +182,7 @@ SUBROUTINE baroclinic_wave_test(deep,moist,pertt,X,lon,lat,p,z,zcoords,u,v,t,the inttau2 = constC * z * exp(- scaledZ**2) ! radius ratio - if (deep .eq. 0) then + if (deep) then rratio = 1.d0 else rratio = (z + aref) / aref; @@ -198,7 +198,7 @@ SUBROUTINE baroclinic_wave_test(deep,moist,pertt,X,lon,lat,p,z,zcoords,u,v,t,the !----------------------------------------------------- inttermU = (rratio * cos(lat))**(K - 1.d0) - (rratio * cos(lat))**(K + 1.d0) bigU = g / aref * K * inttau2 * inttermU * t - if (deep .eq. 0) then + if (deep) then rcoslat = aref * cos(lat) else rcoslat = (z + aref) * cos(lat) @@ -270,7 +270,7 @@ END SUBROUTINE baroclinic_wave_test !----------------------------------------------------------------------- SUBROUTINE evaluate_pressure_temperature(deep, X, lon, lat, z, p, t) - INTEGER, INTENT(IN) :: deep ! Deep (1) or Shallow (0) test case + LOGICAL, INTENT(IN) :: deep ! Deep (1) or Shallow (0) test case REAL(8), INTENT(IN) :: & X, & ! Earth scaling ratio @@ -315,7 +315,7 @@ SUBROUTINE evaluate_pressure_temperature(deep, X, lon, lat, z, p, t) !-------------------------------------------- ! radius ratio !-------------------------------------------- - if (deep .eq. 0) then + if (deep) then rratio = 1.d0 else rratio = (z + aref) / aref; @@ -344,7 +344,7 @@ END SUBROUTINE evaluate_pressure_temperature !----------------------------------------------------------------------- SUBROUTINE evaluate_z_temperature(deep, X, lon, lat, p, z, t) - INTEGER, INTENT(IN) :: deep ! Deep (1) or Shallow (0) test case + LOGICAL, INTENT(IN) :: deep ! Deep (1) or Shallow (0) test case REAL(8), INTENT(IN) :: & X, & ! Earth scaling ratio diff --git a/components/homme/src/test_src/dcmip2016-physics-z.F90 b/components/homme/src/test_src/dcmip2016-physics-z.F90 index fcde94d6ccaa..8962ffe3a79d 100644 --- a/components/homme/src/test_src/dcmip2016-physics-z.F90 +++ b/components/homme/src/test_src/dcmip2016-physics-z.F90 @@ -54,7 +54,7 @@ SUBROUTINE DCMIP2016_PHYSICS(test, u, v, p, theta, qv, qc, qr, rho, & dt, z, zi, lat, nz, precl, pbl_type, prec_type) - use physical_constants, only: g, Rgas, Cp,Rwater_vapor, rearth0, omega0,dd_pi + use physical_constants, only: g=>gravit, Rgas, Cp,Rwater_vapor, rearth0, omega0,dd_pi IMPLICIT NONE diff --git a/components/homme/src/test_src/dcmip2016-supercell.F90 b/components/homme/src/test_src/dcmip2016-supercell.F90 index afd1c2a39140..11fcbf19027a 100644 --- a/components/homme/src/test_src/dcmip2016-supercell.F90 +++ b/components/homme/src/test_src/dcmip2016-supercell.F90 @@ -38,7 +38,7 @@ MODULE supercell ! !======================================================================= - use physical_constants, only: g,p0,kappa,Rgas,Cp,Rwater_vapor,rearth,omega0, dd_pi + use physical_constants, only: g=>gravit,p0,kappa,Rgas,Cp,Rwater_vapor,rearth,omega0, dd_pi IMPLICIT NONE diff --git a/components/homme/src/test_src/dry_planar.F90 b/components/homme/src/test_src/dry_planar.F90 index 3fa5df80105f..b3150ca3395d 100644 --- a/components/homme/src/test_src/dry_planar.F90 +++ b/components/homme/src/test_src/dry_planar.F90 @@ -104,7 +104,7 @@ subroutine general_gravity_wave_init(elem,hybrid,hvcoord,nets,nete,d,f) call gravity_wave(x,y,p,z,zcoords,use_eta,hyam,hybm,d,u,v,w,T,T_mean,phis,ps,rho,rho_mean,q(1)) dp = pressure_thickness(ps,k,hvcoord) call set_state(u,v,w,T,ps,phis,p,dp,zm(k),g, i,j,k,elem(ie),1,nt) - call set_tracers(q,qsize, dp,i,j,k,y,x,elem(ie)) + call set_tracers(q,qsize, dp,zm(k),i,j,k,y,x,elem(ie)) enddo; enddo; enddo; do k=1,nlevp; do j=1,np; do i=1,np call get_xycoordinates(x,y,hyai,hybi, i,j,k,elem(ie),hvcoord) diff --git a/components/homme/src/test_src/held_suarez_mod.F90 b/components/homme/src/test_src/held_suarez_mod.F90 index 7bcd69353834..2ae3ac1603ce 100644 --- a/components/homme/src/test_src/held_suarez_mod.F90 +++ b/components/homme/src/test_src/held_suarez_mod.F90 @@ -13,7 +13,7 @@ module held_suarez_mod use hybrid_mod, only: hybrid_t use hybvcoord_mod, only: hvcoord_t use kinds, only: real_kind, iulog - use physical_constants, only: p0, kappa,g, dd_pi, Rgas + use physical_constants, only: p0, kappa,g=>gravit, dd_pi, Rgas use physics_mod, only: prim_condense use time_mod, only: secpday #ifndef HOMME_WITHOUT_PIOLIBRARY diff --git a/components/homme/src/theta-l/share/element_ops.F90 b/components/homme/src/theta-l/share/element_ops.F90 index fb4d19cec6cb..d1bce092cf8e 100644 --- a/components/homme/src/theta-l/share/element_ops.F90 +++ b/components/homme/src/theta-l/share/element_ops.F90 @@ -59,8 +59,9 @@ module element_ops use kinds, only: real_kind, iulog use perf_mod, only: t_startf, t_stopf, t_barrierf, t_adj_detailf ! _EXTERNAL use parallel_mod, only: abortmp - use physical_constants, only : p0, Cp, Rgas, Rwater_vapor, Cpwater_vapor, kappa, g, dd_pi, TREF + use physical_constants, only : p0, Cp, Rgas, Rwater_vapor, Cpwater_vapor, kappa, rearth, gravit, dd_pi, TREF use control_mod, only: use_moisture, theta_hydrostatic_mode, hv_ref_profiles + use deep_atm_mod, only: quasi_hydrostatic_terms use eos, only: pnh_and_exner_from_eos, phi_from_eos implicit none private @@ -124,7 +125,7 @@ recursive subroutine get_field(elem,name,field,hvcoord,nt,ntQ) if(theta_hydrostatic_mode) then call get_field(elem,'omega',omega,hvcoord,nt,ntQ) call get_field(elem,'rho' ,rho ,hvcoord,nt,ntQ) - field = -omega/(rho *g) + field = -omega/(rho *gravit) else field =elem%state%w_i(:,:,1:nlev,nt) @@ -156,10 +157,11 @@ subroutine get_field_i(elem,name,field,hvcoord,nt) endif case('geo_i'); if(theta_hydrostatic_mode) then - call phi_from_eos(hvcoord,elem%state%phis,elem%state%vtheta_dp(:,:,:,nt),elem%state%dp3d(:,:,:,nt),field) + call phi_from_eos(hvcoord,elem%state%phis,elem%state%ps_v(:,:,nt),elem%state%vtheta_dp(:,:,:,nt),elem%state%dp3d(:,:,:,nt),field) else field = elem%state%phinh_i(:,:,1:nlevp,nt) endif + case('mu_i'); call get_dpnh_dp_i(elem,field,hvcoord,nt) case('pnh_i'); @@ -364,7 +366,7 @@ subroutine get_phi(elem,phi,phi_i,hvcoord,nt) if(theta_hydrostatic_mode) then dp=elem%state%dp3d(:,:,:,nt) - call phi_from_eos(hvcoord,elem%state%phis,elem%state%vtheta_dp(:,:,:,nt),dp,phi_i) + call phi_from_eos(hvcoord,elem%state%phis,elem%state%ps_v(:,:,nt),elem%state%vtheta_dp(:,:,:,nt),dp,phi_i) else phi_i = elem%state%phinh_i(:,:,:,nt) endif @@ -506,6 +508,7 @@ end subroutine set_thermostate !_____________________________________________________________________ subroutine set_state(u,v,w,T,ps,phis,p,dp,zm,g,i,j,k,elem,n0,n1) + use deep_atm_mod, only: r_hat_from_phi, phi_from_z ! ! set state variables at node(i,j,k) at layer midpoints ! used by idealized tests for dry initial conditions @@ -514,19 +517,22 @@ subroutine set_state(u,v,w,T,ps,phis,p,dp,zm,g,i,j,k,elem,n0,n1) real(real_kind), intent(in) :: u,v,w,T,ps,phis,p,dp,zm,g integer, intent(in) :: i,j,k,n0,n1 type(element_t), intent(inout) :: elem + real(real_kind) :: r_hat, pdensity + r_hat = r_hat_from_phi(phi_from_z(zm)) + pdensity = dp ! set prognostic state variables at level midpoints elem%state%v (i,j,1,k,n0:n1) = u elem%state%v (i,j,2,k,n0:n1) = v - elem%state%dp3d(i,j,k, n0:n1) = dp + elem%state%dp3d(i,j,k, n0:n1) = pdensity elem%state%ps_v(i,j, n0:n1) = ps elem%state%phis(i,j) = phis - elem%state%vtheta_dp(i,j,k,n0:n1)=T*dp*((p/p0)**(-kappa)) + elem%state%vtheta_dp(i,j,k,n0:n1)=T*pdensity*((p/p0)**(-kappa)) end subroutine set_state - subroutine set_state_i(u,v,w,T,ps,phis,p,zm,g,i,j,k,elem,n0,n1) + use deep_atm_mod, only: phi_from_z ! ! set state variables at node(i,j,k) at layer interfaces ! used by idealized tests for dry initial conditions @@ -538,12 +544,13 @@ subroutine set_state_i(u,v,w,T,ps,phis,p,zm,g,i,j,k,elem,n0,n1) ! set prognostic state variables at level midpoints elem%state%w_i (i,j,k, n0:n1) = w - elem%state%phinh_i(i,j,k, n0:n1) = g*zm + elem%state%phinh_i(i,j,k, n0:n1) = phi_from_z(zm) end subroutine set_state_i !_____________________________________________________________________ subroutine set_elem_state(u,v,w,w_i,T,ps,phis,p,dp,zm,zi,g,elem,n0,n1,ntQ) + use deep_atm_mod, only: r_hat_from_phi, phi_from_z ! ! set element state variables ! works for both dry and moist initial conditions @@ -557,27 +564,33 @@ subroutine set_elem_state(u,v,w,w_i,T,ps,phis,p,dp,zm,zi,g,elem,n0,n1,ntQ) type(element_t), intent(inout) :: elem integer :: n real(real_kind), dimension(np,np,nlev) :: Rstar + real(real_kind), dimension(np,np,nlev) :: r_hat, pdensity ! get cp and kappa for dry or moist cases call get_R_star(Rstar,elem%state%Q(:,:,:,1)) + + r_hat = r_hat_from_phi(phi_from_z(zm, nlev), nlev) + pdensity = dp + do n=n0,n1 ! set prognostic state variables at level midpoints elem%state%v (:,:,1,:,n) = u elem%state%v (:,:,2,:,n) = v - elem%state%dp3d(:,:,:, n) = dp + elem%state%dp3d(:,:,:, n) = pdensity elem%state%ps_v(:,:, n) = ps elem%state%phis(:,:) = phis - elem%state%vtheta_dp(:,:,:,n) = (Rstar/Rgas)*T*dp*((p/p0)**(-kappa)) + elem%state%vtheta_dp(:,:,:,n) = (Rstar/Rgas)*T*pdensity*((p/p0)**(-kappa)) !/r_hat**2 elem%state%w_i (:,:,:, n) = w_i - elem%state%phinh_i(:,:,:, n) = g*zi + elem%state%phinh_i(:,:,:, n) = phi_from_z(zi, nlevp) end do end subroutine set_elem_state !_____________________________________________________________________ subroutine get_state(u,v,w,T,pnh,dp,ps,rho,zm,zi,g,elem,hvcoord,nt,ntQ) + use deep_atm_mod, only: z_from_phi ! get state variables at layer midpoints ! used by idealized tests to compute idealized physics forcing terms ! currently all forcing is done on u,v and T/theta - no forcing @@ -629,10 +642,10 @@ subroutine get_state(u,v,w,T,pnh,dp,ps,rho,zm,zi,g,elem,hvcoord,nt,ntQ) endif do k=1,nlev - zm(:,:,k) = (phi_i(:,:,k)+phi_i(:,:,k+1))/(2*g) + zm(:,:,k) = (z_from_phi(phi_i(:,:,k))+z_from_phi(phi_i(:,:,k+1)))/(2) end do do k=1,nlevp - zi(:,:,k) = phi_i(:,:,k)/g + zi(:,:,k) = z_from_phi(phi_i(:,:,k)) end do end subroutine get_state @@ -695,7 +708,8 @@ subroutine set_forcing_rayleigh_friction(elem, zm, zi, ztop, zc, tau, u0,v0, n) !_____________________________________________________________________ subroutine tests_finalize(elem,hvcoord,ie) - + use deep_atm_mod, only: r_hat_from_phi, z_from_phi, phi_from_z, g_from_phi + USE, INTRINSIC :: IEEE_ARITHMETIC, ONLY: IEEE_IS_FINITE ! Now that all variables have been initialized, set phi to be in hydrostatic balance implicit none @@ -704,26 +718,33 @@ subroutine tests_finalize(elem,hvcoord,ie) type(element_t), intent(inout):: elem integer, optional, intent(in) :: ie ! optional element index, to save initial state - integer :: k,tl + integer :: k,tl, nind real(real_kind), dimension(np,np,nlev) :: pi - - real(real_kind), dimension(np,np,nlev) :: pnh,exner - real(real_kind), dimension(np,np,nlevp) :: dpnh_dp_i,phi_i + real(real_kind), dimension(np,np,nlev) :: pnh,exner,dphi_tmp, dp + real(real_kind), dimension(np,np,nlevp) :: dpnh_dp_i,phi_i, phi_i_tmp + integer :: niter = 0 + logical, parameter :: do_rootfinding = .false., pnh_instead_of_mu=.false. tl=1 - - call phi_from_eos(hvcoord,elem%state%phis,elem%state%vtheta_dp(:,:,:,tl),& - elem%state%dp3d(:,:,:,tl),elem%state%phinh_i(:,:,:,tl)) + + !call phi_from_eos(hvcoord,elem%state%phis,elem%state%vtheta_dp(:,:,:,tl),& + ! elem%state%dp3d(:,:,:,tl),elem%state%phinh_i(:,:,:,tl)) ! Disable the following check in CUDA bfb builds, ! since the calls to pow are inexact #if !(defined(HOMMEXX_BFB_TESTING) && defined(HOMMEXX_ENABLE_GPU)) - ! verify discrete hydrostatic balance + + + ! verify discrete hydrostatic balance: + ! ===================================== call pnh_and_exner_from_eos(hvcoord,elem%state%vtheta_dp(:,:,:,tl),& elem%state%dp3d(:,:,:,tl),elem%state%phinh_i(:,:,:,tl),pnh,exner,dpnh_dp_i) + + ! ===================== + do k=1,nlev pi(:,:,k) = hvcoord%hyam(k)*hvcoord%ps0 + hvcoord%hybm(k)*elem%state%ps_v(:,:,tl) - if (maxval(abs(1-dpnh_dp_i(:,:,k))) > 1e-10) then + if (maxval(abs(1-dpnh_dp_i(:,:,k))) > 1e-9) then write(iulog,*)'WARNING: hydrostatic inverse FAILED!' write(iulog,*)k,minval(dpnh_dp_i(:,:,k)),maxval(dpnh_dp_i(:,:,k)) write(iulog,*) 'pnh',pi(1,1,k),pnh(1,1,k) @@ -775,7 +796,7 @@ subroutine initialize_reference_states(hvcoord, phis, & ! compute phi_ref temp = theta_ref*dp_ref - call phi_from_eos(hvcoord, phis, temp, dp_ref, phi_ref) + call phi_from_eos(hvcoord, phis, ps_ref, temp, dp_ref, phi_ref) ! keep profiles, based on the value of hv_ref_profiles if (hv_ref_profiles == 0) then @@ -820,7 +841,7 @@ subroutine set_theta_ref(hvcoord,dp,theta_ref,linear_profile) ! reference T = 288K. reference lapse rate = 6.5K/km = .0065 K/m ! Tref = T0+T1*exner ! Thetaref = T0/exner + T1 - T1 = tref_lapse_rate*TREF*Cp/g ! = 191 + T1 = tref_lapse_rate*TREF*Cp/gravit ! = 191 WARNING: NOT ADAPTED TO DEEP ATMOSPHERE T0 = TREF-T1 ! = 97 p_i(:,:,1) = hvcoord%hyai(1)*hvcoord%ps0 diff --git a/components/homme/src/theta-l/share/eos.F90 b/components/homme/src/theta-l/share/eos.F90 index faad7ba744e2..45b6f6dbe677 100644 --- a/components/homme/src/theta-l/share/eos.F90 +++ b/components/homme/src/theta-l/share/eos.F90 @@ -22,7 +22,7 @@ module eos use hybvcoord_mod, only: hvcoord_t use kinds, only: real_kind use parallel_mod, only: abortmp - use physical_constants, only : p0, kappa, g, Rgas + use physical_constants, only : p0, kappa, gravit, Rgas, rearth use control_mod, only: theta_hydrostatic_mode #ifdef HOMMEXX_BFB_TESTING use bfb_mod, only: bfb_pow @@ -59,6 +59,7 @@ subroutine pnh_and_exner_from_eos(hvcoord,vtheta_dp,dp3d,phi_i,pnh,exner,& real (kind=real_kind), intent(out), optional :: pnh_i_out(np,np,nlevp) ! pnh on interfaces character(len=*), intent(in), optional :: caller ! name for error + ! local real (kind=real_kind) :: dphi(np,np,nlev) integer :: k @@ -66,19 +67,21 @@ subroutine pnh_and_exner_from_eos(hvcoord,vtheta_dp,dp3d,phi_i,pnh,exner,& do k=1,nlev dphi(:,:,k)=phi_i(:,:,k+1)-phi_i(:,:,k) enddo + !write(*,*) "dphi", dphi(1,1,:) if (present(caller)) then - call pnh_and_exner_from_eos2(hvcoord,vtheta_dp,dp3d,dphi,pnh,exner,& + call pnh_and_exner_from_eos2(hvcoord,vtheta_dp,dp3d,dphi,phi_i,pnh,exner,& dpnh_dp_i,caller,pnh_i_out) else - call pnh_and_exner_from_eos2(hvcoord,vtheta_dp,dp3d,dphi,pnh,exner,& + call pnh_and_exner_from_eos2(hvcoord,vtheta_dp,dp3d,dphi,phi_i,pnh,exner,& dpnh_dp_i,'not specified',pnh_i_out) endif end subroutine pnh_and_exner_from_eos -subroutine pnh_and_exner_from_eos2(hvcoord,vtheta_dp,dp3d,dphi,pnh,exner,& +subroutine pnh_and_exner_from_eos2(hvcoord,vtheta_dp,dp3d,dphi,phi_i,pnh,exner,& dpnh_dp_i,caller,pnh_i_out) +use deep_atm_mod, only: r_hat_from_phi, g_from_phi implicit none ! ! Use Equation of State to compute exner pressure, nh presure @@ -98,10 +101,11 @@ subroutine pnh_and_exner_from_eos2(hvcoord,vtheta_dp,dp3d,dphi,pnh,exner,& real (kind=real_kind), intent(in) :: vtheta_dp(np,np,nlev) real (kind=real_kind), intent(in) :: dp3d(np,np,nlev) real (kind=real_kind), intent(in) :: dphi(np,np,nlev) + real (kind=real_kind), intent(in) :: phi_i(np,np,nlevp) ! Needed for deep atmosphere real (kind=real_kind), intent(out) :: pnh(np,np,nlev) ! nh nonhyrdo pressure real (kind=real_kind), intent(out) :: dpnh_dp_i(np,np,nlevp) ! d(pnh) / d(pi) real (kind=real_kind), intent(out) :: exner(np,np,nlev) ! exner nh pressure - character(len=*), intent(in) :: caller ! name for error + character(len=*), intent(in), optional :: caller ! name for error real (kind=real_kind), intent(out), optional :: pnh_i_out(np,np,nlevp) ! pnh on interfaces @@ -112,6 +116,7 @@ subroutine pnh_and_exner_from_eos2(hvcoord,vtheta_dp,dp3d,dphi,pnh,exner,& real (kind=real_kind) :: pnh_i(np,np,nlevp) real (kind=real_kind) :: dp3d_i(np,np,nlevp) real (kind=real_kind) :: pi_i(np,np,nlevp) + real (kind=real_kind) :: r_hat(np,np) integer :: i,j,k,k2 logical :: ierr @@ -173,7 +178,9 @@ subroutine pnh_and_exner_from_eos2(hvcoord,vtheta_dp,dp3d,dphi,pnh,exner,& ! non-hydrostatic EOS !============================================================== do k=1,nlev - p_over_exner(:,:,k) = Rgas*vtheta_dp(:,:,k)/(-dphi(:,:,k)) + ! DEEP ATMOSPHERE MODIFICATION + r_hat = r_hat_from_phi((phi_i(:, :, k) + phi_i(:, :, k+1))/2) ! DA_CHANGE + p_over_exner(:,:,k) = Rgas*vtheta_dp(:,:,k)/(-dphi(:,:,k))/r_hat**2 #ifndef HOMMEXX_BFB_TESTING pnh(:,:,k) = p0 * (p_over_exner(:,:,k)/p0)**(1/(1-kappa)) #else @@ -184,12 +191,15 @@ subroutine pnh_and_exner_from_eos2(hvcoord,vtheta_dp,dp3d,dphi,pnh,exner,& !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! boundary terms !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - pnh_i(:,:,1) = hvcoord%hyai(1)*hvcoord%ps0 ! hydrostatic ptop + r_hat = r_hat_from_phi(phi_i(:, :, 1)) + !pnh_i(:,:,1) = 0.05 * hvcoord%hyai(1)*hvcoord%ps0 ! hydrostatic ptop + !pnh_i(:,:,1) = pnh(:,:,1)-dp3d(:,:,1)/2 ! hydrostatic ptop ! surface boundary condition pnh_i determined by w equation to enforce ! w b.c. This is computed in the RHS calculation. Here, we use ! an approximation (hydrostatic) so that dpnh/dpi = 1 ! DO NOT CHANGE this approximation. it is required by ! compute_andor_apply_rhs() + pnh_i(:,:,1) = pnh(:,:,1) - dp3d(:, :,1)/(2.0 * r_hat**2) !* (1 + u_top/g_from_phi(phi_i(:, :, 1))) pnh_i(:,:,nlevp) = pnh(:,:,nlev) + dp3d(:,:,nlev)/2 @@ -200,14 +210,14 @@ subroutine pnh_and_exner_from_eos2(hvcoord,vtheta_dp,dp3d,dphi,pnh,exner,& do k=2,nlev dp3d_i(:,:,k)=(dp3d(:,:,k)+dp3d(:,:,k-1))/2 end do - - dpnh_dp_i(:,:,1) = 2*(pnh(:,:,1)-pnh_i(:,:,1))/dp3d_i(:,:,1) - dpnh_dp_i(:,:,nlevp) = 2*(pnh_i(:,:,nlevp)-pnh(:,:,nlev))/dp3d_i(:,:,nlevp) + r_hat = r_hat_from_phi(phi_i(:, :, 1)) ! DA_CHANGE + dpnh_dp_i(:,:,1) = 2*(pnh(:,:,1)-pnh_i(:,:,1))/dp3d_i(:,:,1) * r_hat**2 + r_hat = r_hat_from_phi(phi_i(:, :, nlevp)) ! DA_CHANGE + dpnh_dp_i(:,:,nlevp) = 2*(pnh_i(:,:,nlevp)-pnh(:,:,nlev))/dp3d_i(:,:,nlevp) * r_hat**2 do k=2,nlev - dpnh_dp_i(:,:,k) = (pnh(:,:,k)-pnh(:,:,k-1))/dp3d_i(:,:,k) + r_hat = r_hat_from_phi(phi_i(:, :, k)) !DA_CHANGE + dpnh_dp_i(:,:,k) = (pnh(:,:,k) -pnh(:,:,k-1) )/dp3d_i(:,:,k) * r_hat**2 end do - - if (present(pnh_i_out)) then ! boundary values already computed. interpolate interior ! use linear interpolation in hydrostatic pressure coordinate @@ -225,7 +235,8 @@ subroutine pnh_and_exner_from_eos2(hvcoord,vtheta_dp,dp3d,dphi,pnh,exner,& !_____________________________________________________________________ - subroutine phi_from_eos(hvcoord,phis,vtheta_dp,dp,phi_i) + subroutine phi_from_eos(hvcoord,phis,ps,vtheta_dp,dp,phi_i) + use deep_atm_mod, only: r_hat_from_phi ! ! Use Equation of State to compute geopotential ! @@ -249,12 +260,14 @@ subroutine phi_from_eos(hvcoord,phis,vtheta_dp,dp,phi_i) type (hvcoord_t), intent(in) :: hvcoord ! hybrid vertical coordinate struct real (kind=real_kind), intent(in) :: vtheta_dp(np,np,nlev) real (kind=real_kind), intent(in) :: dp(np,np,nlev) - real (kind=real_kind), intent(in) :: phis(np,np) + real (kind=real_kind), intent(in) :: phis(np,np), ps(np,np) real (kind=real_kind), intent(out) :: phi_i(np,np,nlevp) ! local + real (kind=real_kind) :: dphi_guess(np,np) ! pressure at cell centers real (kind=real_kind) :: p(np,np,nlev) ! pressure at cell centers real (kind=real_kind) :: p_i(np,np,nlevp) ! pressure on interfaces + real (kind=real_kind) :: r_hat(np,np) ! pressure at cell centers integer :: k @@ -284,9 +297,9 @@ subroutine phi_from_eos(hvcoord,phis,vtheta_dp,dp,phi_i) endif #endif ! compute pressure on interfaces - p_i(:,:,1)=hvcoord%hyai(1)*hvcoord%ps0 - do k=1,nlev - p_i(:,:,k+1)=p_i(:,:,k) + dp(:,:,k) + p_i(:,:,nlevp)=ps + do k=nlev,1,-1 + p_i(:,:,k)=p_i(:,:,k+1) - dp(:,:,k) enddo do k=1,nlev p(:,:,k) = (p_i(:,:,k+1)+p_i(:,:,k))/2 @@ -297,7 +310,10 @@ subroutine phi_from_eos(hvcoord,phis,vtheta_dp,dp,phi_i) #ifdef HOMMEXX_BFB_TESTING phi_i(:,:,k) = phi_i(:,:,k+1)+ (Rgas*vtheta_dp(:,:,k)*bfb_pow(p(:,:,k)/p0,(kappa-1)))/p0 #else - phi_i(:,:,k) = phi_i(:,:,k+1)+(Rgas*vtheta_dp(:,:,k)*(p(:,:,k)/p0)**(kappa-1))/p0 + r_hat = r_hat_from_phi(phi_i(:,:,k+1)) + dphi_guess = (Rgas*vtheta_dp(:,:,k)*(p(:,:,k)/p0)**(kappa-1))/p0 + r_hat = r_hat_from_phi(phi_i(:,:,k+1) + dphi_guess/r_hat**2) + phi_i(:,:,k) = phi_i(:,:,k+1)+dphi_guess/r_hat**2 #endif enddo end subroutine diff --git a/components/homme/src/theta-l/share/imex_mod.F90 b/components/homme/src/theta-l/share/imex_mod.F90 index a43ca6fb6aea..b6898be7254e 100644 --- a/components/homme/src/theta-l/share/imex_mod.F90 +++ b/components/homme/src/theta-l/share/imex_mod.F90 @@ -10,7 +10,7 @@ module imex_mod use element_mod, only: element_t use derivative_mod, only: derivative_t use time_mod, only: timelevel_t, timelevel_qdp - use physical_constants, only: g, kappa + use physical_constants, only: gravit, kappa, rearth use eos, only: pnh_and_exner_from_eos, pnh_and_exner_from_eos2, phi_from_eos use element_state, only: max_itercnt, max_deltaerr, max_reserr use control_mod, only: theta_hydrostatic_mode, qsplit @@ -59,6 +59,7 @@ end subroutine compute_gwphis subroutine compute_stage_value_dirk(nm1,alphadt_nm1,n0,alphadt_n0,np1,dt2,elem,hvcoord,hybrid,& deriv,nets,nete,itercount,itererr,verbosity_in) + use deep_atm_mod, only: r_hat_from_phi, g_from_phi, z_from_phi !=================================================================================== ! this subroutine solves a stage value equation for a DIRK method which takes the form ! @@ -116,6 +117,7 @@ subroutine compute_stage_value_dirk(nm1,alphadt_nm1,n0,alphadt_n0,np1,dt2,elem,h real (kind=real_kind) :: Fn(np,np,nlev),x(np,np,nlev) real (kind=real_kind) :: gwh_i(np,np,nlevp) ! w hydrostatic real (kind=real_kind) :: v_i(np,np,2,nlevp) ! w hydrostatic + real (kind=real_kind) :: recip_r_hat(np,np,nlev) real (kind=real_kind) :: Jac2D(np,np,nlev) , Jac2L(np,np,nlev-1) real (kind=real_kind) :: Jac2U(np,np,nlev-1) @@ -143,7 +145,7 @@ subroutine compute_stage_value_dirk(nm1,alphadt_nm1,n0,alphadt_n0,np1,dt2,elem,h itercount=0 itererr=0 do ie=nets,nete - call phi_from_eos(hvcoord,elem(ie)%state%phis,elem(ie)%state%vtheta_dp(:,:,:,np1),& + call phi_from_eos(hvcoord,elem(ie)%state%phis, elem(ie)%state%ps_v(:,:,np1),elem(ie)%state%vtheta_dp(:,:,:,np1),& elem(ie)%state%dp3d(:,:,:,np1),elem(ie)%state%phinh_i(:,:,:,np1)) elem(ie)%state%w_i(:,:,:,np1)=0 enddo @@ -179,13 +181,13 @@ subroutine compute_stage_value_dirk(nm1,alphadt_nm1,n0,alphadt_n0,np1,dt2,elem,h call pnh_and_exner_from_eos(hvcoord,elem(ie)%state%vtheta_dp(:,:,:,nt), & elem(ie)%state%dp3d(:,:,:,nt),elem(ie)%state%phinh_i(:,:,:,nt),pnh, & exner,dpnh_dp_i,caller='dirkn0') - w_n0(:,:,1:nlev) = w_n0(:,:,1:nlev) + & - dt3*g*(dpnh_dp_i(:,:,1:nlev)-1) + recip_r_hat = 1.0_real_kind/r_hat_from_phi(elem(ie)%state%phinh_i(:,:,1:nlev,nt), nlev)! DA_CHANGE + w_n0(:,:,1:nlev) = w_n0(:,:,1:nlev) + dt3*(g_from_phi(elem(ie)%state%phinh_i(:,:,1:nlev,nt), nlev))*(dpnh_dp_i(:,:,1:nlev)-1) ! DA_CHANGE call compute_gwphis(gwh_i,elem(ie)%state%dp3d(:,:,:,nt),elem(ie)%state%v(:,:,:,:,nt),& elem(ie)%derived%gradphis,hvcoord) phi_n0(:,:,1:nlev) = phi_n0(:,:,1:nlev) + & - dt3*g*elem(ie)%state%w_i(:,:,1:nlev,nt) - dt3*gwh_i(:,:,1:nlev) + dt3*(g_from_phi(elem(ie)%state%phinh_i(:,:,1:nlev,nt),nlev))*elem(ie)%state%w_i(:,:,1:nlev,nt) - dt3*gwh_i(:,:,1:nlev) !DA_CHANGE end if @@ -196,12 +198,12 @@ subroutine compute_stage_value_dirk(nm1,alphadt_nm1,n0,alphadt_n0,np1,dt2,elem,h elem(ie)%state%dp3d(:,:,:,nt),elem(ie)%state%phinh_i(:,:,:,nt),pnh, & exner,dpnh_dp_i,caller='dirknm1') w_n0(:,:,1:nlev) = w_n0(:,:,1:nlev) + & - dt3*g*(dpnh_dp_i(:,:,1:nlev)-1) + dt3*(g_from_phi(elem(ie)%state%phinh_i(:,:,1:nlev,nt),nlev))*(dpnh_dp_i(:,:,1:nlev)-1) ! DA_CHANGE call compute_gwphis(gwh_i,elem(ie)%state%dp3d(:,:,:,nt),elem(ie)%state%v(:,:,:,:,nt),& elem(ie)%derived%gradphis,hvcoord) phi_n0(:,:,1:nlev) = phi_n0(:,:,1:nlev) + & - dt3*g*elem(ie)%state%w_i(:,:,1:nlev,nt) - dt3*gwh_i(:,:,1:nlev) + dt3*g_from_phi(elem(ie)%state%phinh_i(:,:,1:nlev,nt),nlev)*elem(ie)%state%w_i(:,:,1:nlev,nt) - dt3*gwh_i(:,:,1:nlev) ! DA_CHANGE end if ! add just the wphis(np1) term to the RHS @@ -225,9 +227,9 @@ subroutine compute_stage_value_dirk(nm1,alphadt_nm1,n0,alphadt_n0,np1,dt2,elem,h #endif #if 1 ! use hydrostatic for initial guess - call phi_from_eos(hvcoord,elem(ie)%state%phis,elem(ie)%state%vtheta_dp(:,:,:,np1),& + call phi_from_eos(hvcoord,elem(ie)%state%phis,elem(ie)%state%ps_v(:,:,np1),elem(ie)%state%vtheta_dp(:,:,:,np1),& elem(ie)%state%dp3d(:,:,:,np1),phi_np1) - w_np1(:,:,1:nlev) = (phi_np1(:,:,1:nlev) - phi_n0(:,:,1:nlev) )/(dt2*g) + w_np1(:,:,1:nlev) = (z_from_phi(phi_np1(:,:,1:nlev), nlev) - z_from_phi(phi_n0(:,:,1:nlev), nlev ))/(dt2) ! DA_CHANGE #endif ! initial residual @@ -240,12 +242,12 @@ subroutine compute_stage_value_dirk(nm1,alphadt_nm1,n0,alphadt_n0,np1,dt2,elem,h do k=1,nlev do j=1,np do i=1,np - if ( dphi(i,j,k) > -g) then + if ( dphi(i,j,k) > -gravit) then ! DA_CHANGE WARNING: not adapted to Deep Atmosphere if (verbosity > 0) then write(iulog,*) 'WARNING:IMEX limiting initial guess. ie,i,j,k=',ie,i,j,k write(iulog,*) 'dphi(i,j,k)= ',dphi(i,j,k) endif - dphi(i,j,k)=-g + dphi(i,j,k)=-gravit ! DA_CHANGE nsafe=1 endif enddo @@ -256,20 +258,24 @@ subroutine compute_stage_value_dirk(nm1,alphadt_nm1,n0,alphadt_n0,np1,dt2,elem,h do k=nlev,1,-1 ! scan phi_np1(:,:,k) = phi_np1(:,:,k+1)-dphi(:,:,k) enddo - w_np1(:,:,1:nlev) = (phi_np1(:,:,1:nlev) - phi_n0(:,:,1:nlev) )/(dt2*g) + w_np1(:,:,1:nlev) = (z_from_phi(phi_np1(:,:,1:nlev), nlev) - z_from_phi(phi_n0(:,:,1:nlev), nlev) )/(dt2) !DA_CHANGE endif call pnh_and_exner_from_eos2(hvcoord,elem(ie)%state%vtheta_dp(:,:,:,np1),elem(ie)%state%dp3d(:,:,:,np1),& - dphi,pnh,exner,dpnh_dp_i,'dirk1') + dphi,phi_np1,pnh,exner,dpnh_dp_i,'dirk1') Fn(:,:,1:nlev) = w_np1(:,:,1:nlev) - & - (w_n0(:,:,1:nlev) + g*dt2 * (dpnh_dp_i(:,:,1:nlev)-1)) + (w_n0(:,:,1:nlev) + g_from_phi(elem(ie)%state%phinh_i(:,:,1:nlev,np1),nlev)*dt2 * (dpnh_dp_i(:,:,1:nlev)-1)) itercount=0 do while (itercount < maxiter) + do k=nlev,1,-1 ! scan + phi_np1(:,:,k) = phi_np1(:,:,k+1)-dphi(:,:,k) + enddo + ! numerical J: - !call get_dirk_jacobian(JacL,JacD,JacU,dt2,elem(ie)%state%dp3d(:,:,:,np1),dphi,pnh,0,1d-4,hvcoord,dpnh_dp_i,vtheta_dp) + !call get_dirk_jacobian(JacL,JacD,JacU,dt2,elem(ie)%state%dp3d(:,:,:,np1),dphi, phi_np1,pnh,0,1d-4,hvcoord,dpnh_dp_i,elem(ie)%state%vtheta_dp(:,:,:,np1)) ! analytic J: - call get_dirk_jacobian(JacL,JacD,JacU,dt2,elem(ie)%state%dp3d(:,:,:,np1),dphi,pnh,1) + call get_dirk_jacobian(JacL,JacD,JacU,dt2,elem(ie)%state%dp3d(:,:,:,np1),dphi,phi_np1,pnh,1) x(:,:,1:nlev) = -Fn(:,:,1:nlev) @@ -277,10 +283,11 @@ subroutine compute_stage_value_dirk(nm1,alphadt_nm1,n0,alphadt_n0,np1,dt2,elem,h do k = 1,nlev-1 dphi(:,:,k) = dphi_n0(:,:,k) + & - dt2*g*((w_np1(:,:,k+1) - w_np1(:,:,k)) + & + dt2*g_from_phi((elem(ie)%state%phinh_i(:,:,k,nt) + elem(ie)%state%phinh_i(:,:,k+1,nt))/2)*((w_np1(:,:,k+1) - w_np1(:,:,k)) + & !DA_CHANGE (x(:,:,k+1) - x(:,:,k))) end do - dphi(:,:,nlev) = dphi_n0(:,:,nlev) - dt2*g*(w_np1(:,:,nlev) + x(:,:,nlev)) + dphi(:,:,nlev) = dphi_n0(:,:,nlev) - dt2*g_from_phi((elem(ie)%state%phinh_i(:,:,nlev,nt) + & + elem(ie)%state%phinh_i(:,:,nlev+1,nt))/2)*(w_np1(:,:,nlev) + x(:,:,nlev)) ! DA_CHANGE alphas = 1 if (any(dphi(:,:,1:nlev) >= 0)) then @@ -299,7 +306,9 @@ subroutine compute_stage_value_dirk(nm1,alphadt_nm1,n0,alphadt_n0,np1,dt2,elem,h dw = -w_np1(i,j,k) end if if (dx /= 0) then - alpha_k = -(dphi_n0(i,j,k) + dt2*g*dw)/(dt2*g*dx) + alpha_k = -(dphi_n0(i,j,k) + dt2*g_from_phi((elem(ie)%state%phinh_i(i,j,k,nt)+ & + elem(ie)%state%phinh_i(i,j,k+1,nt))/2.0)*dw)/(dt2*g_from_phi((elem(ie)%state%phinh_i(i,j,k,nt)+ & + elem(ie)%state%phinh_i(i,j,k+1,nt))/2.0)*dx) ! DA_CHANGE if (alpha_k >= 0) alpha = min(alpha, alpha_k) end if end do @@ -310,25 +319,28 @@ subroutine compute_stage_value_dirk(nm1,alphadt_nm1,n0,alphadt_n0,np1,dt2,elem,h do k = 1,nlev-1 dphi(i,j,k) = dphi_n0(i,j,k) + & - dt2*g*((w_np1(i,j,k+1) - w_np1(i,j,k)) + & - alpha*(x(i,j,k+1) - x(i,j,k))) + dt2*g_from_phi((elem(ie)%state%phinh_i(i,j,k,nt)+& + elem(ie)%state%phinh_i(i,j,k+1,nt))/2.0)*((w_np1(i,j,k+1) - w_np1(i,j,k)) + alpha*(x(i,j,k+1) - x(i,j,k))) !DA_CHANGE end do - dphi(i,j,nlev) = dphi_n0(i,j,nlev) - dt2*g*(w_np1(i,j,nlev) + alpha*x(i,j,nlev)) + dphi(i,j,nlev) = dphi_n0(i,j,nlev) - dt2*g_from_phi((elem(ie)%state%phinh_i(i,j,k,nt)+ elem(ie)%state%phinh_i(i,j,k+1,nt))/2.0)*(w_np1(i,j,nlev) + alpha*x(i,j,nlev)) !DA_CHANGE end do end do end if do k = 1,nlev w_np1(:,:,k) = w_np1(:,:,k) + alphas*x(:,:,k) end do - + do k=nlev,1,-1 ! scan + phi_np1(:,:,k) = phi_np1(:,:,k+1)-dphi(:,:,k) + enddo + call pnh_and_exner_from_eos2(hvcoord,elem(ie)%state%vtheta_dp(:,:,:,np1),& - elem(ie)%state%dp3d(:,:,:,np1),dphi,pnh,exner,dpnh_dp_i,'dirk2') - Fn(:,:,1:nlev) = w_np1(:,:,1:nlev) - (w_n0(:,:,1:nlev) + g*dt2 * (dpnh_dp_i(:,:,1:nlev)-1)) + elem(ie)%state%dp3d(:,:,:,np1),dphi,phi_np1,pnh,exner,dpnh_dp_i,'dirk2') + Fn(:,:,1:nlev) = w_np1(:,:,1:nlev) - (w_n0(:,:,1:nlev) + g_from_phi((elem(ie)%state%phinh_i(:,:,1:nlev,nt)+ & + elem(ie)%state%phinh_i(:,:,2:nlevp,nt))/2.0, nlev)*dt2 * (dpnh_dp_i(:,:,1:nlev)-1)) ! DA_CHANGE ! this is not used in this loop, so move out of loop !reserr=maxval(abs(Fn))/(wmax*abs(dt2)) deltaerr=maxval(abs(x))/wmax - ! update iteration count and error measure itercount=itercount+1 !if (reserr < restol) exit @@ -336,7 +348,7 @@ subroutine compute_stage_value_dirk(nm1,alphadt_nm1,n0,alphadt_n0,np1,dt2,elem,h end do ! end do for the do while loop ! update phi: - phi_np1(:,:,1:nlev) = phi_n0(:,:,1:nlev) + dt2*g*w_np1(:,:,1:nlev) + phi_np1(:,:,1:nlev) = phi_n0(:,:,1:nlev) + dt2*g_from_phi(phi_n0(:,:,1:nlev), nlev)*w_np1(:,:,1:nlev) !DA_CHANGE ! keep track of running max iteraitons and max error (reset after each diagnostics output) @@ -364,8 +376,9 @@ subroutine compute_stage_value_dirk(nm1,alphadt_nm1,n0,alphadt_n0,np1,dt2,elem,h end subroutine compute_stage_value_dirk - subroutine get_dirk_jacobian(JacL,JacD,JacU,dt2,dp3d,dphi,pnh,exact,& + subroutine get_dirk_jacobian(JacL,JacD,JacU,dt2,dp3d,dphi,phi_i,pnh,exact,& epsie,hvcoord,dpnh_dp_i,vtheta_dp) + use deep_atm_mod, only: r_hat_from_phi, g_from_phi !================================================================================ ! compute Jacobian of F(phi) = sum(dphi) +const + (dt*g)^2 *(1-dp/dpi) column wise ! with respect to phi @@ -390,6 +403,7 @@ subroutine get_dirk_jacobian(JacL,JacD,JacU,dt2,dp3d,dphi,pnh,exact,& real (kind=real_kind), intent(in) :: dp3d(np,np,nlev) real (kind=real_kind), intent(inout) :: pnh(np,np,nlev) real (kind=real_kind), intent(in) :: dphi(np,np,nlev) + real (kind=real_kind), intent(in) :: phi_i(np,np,nlevp) real (kind=real_kind), intent(in) :: dt2 real (kind=real_kind), intent(in), optional :: epsie ! epsie is the differencing size in the approx. Jacobian @@ -401,36 +415,44 @@ subroutine get_dirk_jacobian(JacL,JacD,JacU,dt2,dp3d,dphi,pnh,exact,& ! local real (kind=real_kind) :: alpha1(np,np),alpha2(np,np) - real (kind=real_kind) :: e(np,np,nlev),dphi_temp(np,np,nlev),exner(np,np,nlev) + real (kind=real_kind) :: e(np,np,nlev),dphi_temp(np,np,nlev),exner(np,np,nlev),phi_i_tmp(np,np,nlevp) real (kind=real_kind) :: dpnh2(np,np,nlev),dpnh_dp_i_epsie(np,np,nlevp) - real (kind=real_kind) :: ds(np,np,nlev),delta_mu(np,np,nlevp) - real (kind=real_kind) :: a,b(np,np),ck(np,np),ckm1(np,np) + real (kind=real_kind) :: ds(np,np,nlev),delta_mu(np,np,nlevp),r_hat(np,np,nlevp) + real (kind=real_kind) :: a(np,np),b(np,np),ck(np,np),ckm1(np,np) ! integer :: k,l,k2 - if (exact.eq.1) then ! use exact Jacobian + if (exact.eq.1) then ! use exact JacobianA + ! this code will need to change when the equation of state is changed. ! add special cases for k==1 and k==nlev+1 - - a = (dt2*g)**2/(1-kappa) + phi_i_tmp = phi_i k = 1 ! Jacobian row 1 + r_hat(:,:,k) = r_hat_from_phi(phi_i_tmp(:,:,k)) ! DA_CHANGE + a = (dt2*(g_from_phi(phi_i_tmp(:,:,k))))**2/(1-kappa) ! DA_CHANGE b = a/dp3d(:,:,k) ck = pnh(:,:,k)/dphi(:,:,k) - JacU(:,:,k) = 2*b*ck + JacU(:,:,k) = 0 * r_hat(:,:,k)**2 * 2*b*ck ! Neumann upper boundary sets \mu \equiv 1 JacD(:,:,k) = 1 - JacU(:,:,k) ckm1 = ck do k = 2,nlev-1 ! Jacobian row k + r_hat(:,:,k) = r_hat_from_phi(phi_i_tmp(:,:,k)) ! DA_CHANGE + a = (dt2*(g_from_phi(phi_i_tmp(:,:,k))))**2/(1-kappa) ! DA_CHANGE b = 2*a/(dp3d(:,:,k-1) + dp3d(:,:,k)) ck = pnh(:,:,k)/dphi(:,:,k) - JacL(:,:,k-1) = b*ckm1 - JacU(:,:,k ) = b*ck + JacL(:,:,k-1) = r_hat(:,:,k)**2 * b*ckm1 ! TODO: this disagrees with definition of b that comes from Neumann + JacU(:,:,k ) = r_hat(:,:,k)**2 * b*ck JacD(:,:,k ) = 1 - JacL(:,:,k-1) - JacU(:,:,k) ckm1 = ck end do + + !JacL(:,:,1) = 0 k = nlev ! Jacobian row nlev + r_hat(:,:,k) = r_hat_from_phi( phi_i_tmp(:,:,k)) !DA_CHANGE + a = (dt2*(g_from_phi(phi_i_tmp(:,:,k))))**2/(1-kappa) !DA_CHANGE b = 2*a/(dp3d(:,:,k) + dp3d(:,:,k-1)) ck = pnh(:,:,k)/dphi(:,:,k) - JacL(:,:,k-1) = b*ckm1 - JacD(:,:,k ) = 1 - JacL(:,:,k-1) - b*ck + JacL(:,:,k-1) = r_hat(:,:,k)**2 * b*ckm1 + JacD(:,:,k ) = 1 - JacL(:,:,k-1) - r_hat(:,:,k)**2 * b*ck else ! use finite difference approximation to Jacobian with differencing size espie ! compute Jacobian of F(dphi) = phi +const + (dt*g)^2 *(1-dp/dpi) column wise @@ -439,6 +461,7 @@ subroutine get_dirk_jacobian(JacL,JacD,JacU,dt2,dp3d,dphi,pnh,exact,& ! we only form the tridagonal entries and this code can easily be modified to ! accomodate sparse non-tridigonal and dense Jacobians, however, testing only ! the tridiagonal of a Jacobian is probably sufficient for testing purpose + phi_i_tmp = phi_i do k=1,nlev e=0 e(:,:,k)=1 @@ -452,8 +475,13 @@ subroutine get_dirk_jacobian(JacL,JacD,JacU,dt2,dp3d,dphi,pnh,exact,& !dpnh_dp_i_epsie(:,:,:)=1.d0 delta_mu=0 else - call pnh_and_exner_from_eos2(hvcoord,vtheta_dp,dp3d,dphi_temp,pnh,exner,dpnh_dp_i_epsie,'get_dirk_jacobian') - delta_mu(:,:,:)=(g*dt2)**2*(dpnh_dp_i(:,:,:)-dpnh_dp_i_epsie(:,:,:))/epsie + do k2=nlev,1,-1 ! scan + phi_i_tmp(:,:,k2) = phi_i_tmp(:,:,k2+1)-dphi_temp(:,:,k2) + enddo + + call pnh_and_exner_from_eos2(hvcoord,vtheta_dp,dp3d,dphi_temp,phi_i_tmp, pnh,exner,dpnh_dp_i_epsie,'get_dirk_jacobian') + r_hat = r_hat_from_phi(phi_i_tmp, nlevp) !DA_CHANGE + delta_mu(:,:,:)= (g_from_phi(phi_i_tmp,nlevp)*dt2)**2*(dpnh_dp_i(:,:,:)-dpnh_dp_i_epsie(:,:,:))/epsie ! DA_CHANGE end if JacD(:,:,k) = 1 + delta_mu(:,:,k) @@ -510,7 +538,7 @@ subroutine test_imex_jacobian(elem,hybrid,hvcoord,tl,nets,nete) do k=1,nlev dphi(:,:,k)=phi_i(:,:,k+1)-phi_i(:,:,k) enddo - call get_dirk_jacobian(JacL,JacD,JacU,dt,dp3d,dphi,pnh,1) + call get_dirk_jacobian(JacL,JacD,JacU,dt,dp3d,dphi,phi_i,pnh,1) ! compute infinity norm of the initial Jacobian norminfJ0=0.d0 @@ -546,7 +574,7 @@ subroutine test_imex_jacobian(elem,hybrid,hvcoord,tl,nets,nete) do k=1,nlev dphi(:,:,k)=phi_i(:,:,k+1)-phi_i(:,:,k) enddo - call get_dirk_jacobian(Jac2L,Jac2D,Jac2U,dt,dp3d,dphi,pnh,0,& + call get_dirk_jacobian(Jac2L,Jac2D,Jac2U,dt,dp3d,dphi,phi_i,pnh,0,& epsie,hvcoord,dpnh_dp_i,vtheta_dp) if (maxval(abs(JacD(:,:,:)-Jac2D(:,:,:))) > jacerrorvec(j)) then diff --git a/components/homme/src/theta-l/share/model_init_mod.F90 b/components/homme/src/theta-l/share/model_init_mod.F90 index 860508b12529..b5b888b66a11 100644 --- a/components/homme/src/theta-l/share/model_init_mod.F90 +++ b/components/homme/src/theta-l/share/model_init_mod.F90 @@ -25,7 +25,7 @@ module model_init_mod use control_mod, only: qsplit,theta_hydrostatic_mode, hv_ref_profiles, & hv_theta_correction, tom_sponge_start use time_mod, only: timelevel_qdp, timelevel_t - use physical_constants, only: g, TREF, Rgas, kappa + use physical_constants, only: gravit, TREF, Rgas, kappa, rearth use imex_mod, only: test_imex_jacobian use eos, only: phi_from_eos @@ -34,6 +34,7 @@ module model_init_mod contains subroutine model_init2(elem,hybrid,deriv,hvcoord,tl,nets,nete ) + use deep_atm_mod, only: r_hat_from_phi, g_from_phi type(element_t) , intent(inout) :: elem(:) type(hybrid_t) , intent(in) :: hybrid @@ -45,13 +46,17 @@ subroutine model_init2(elem,hybrid,deriv,hvcoord,tl,nets,nete ) ! local variables integer :: ie,t,k real (kind=real_kind) :: gradtemp(np,np,2,nets:nete) + real (kind=real_kind) :: recip_r_hat(np, np) real (kind=real_kind) :: temp(np,np,nlev),ps_ref(np,np) real (kind=real_kind) :: ptop_over_press ! other theta specific model initialization should go here do ie=nets,nete + recip_r_hat = 1.0_real_kind/r_hat_from_phi(elem(ie)%state%phis(:,:)) !DA_CHANGE gradtemp(:,:,:,ie) = gradient_sphere( elem(ie)%state%phis(:,:), deriv, elem(ie)%Dinv) + gradtemp(:,:,1,ie) = recip_r_hat * gradtemp(:,:,1,ie) + gradtemp(:,:,2,ie) = recip_r_hat * gradtemp(:,:,2,ie) enddo call make_C0_vector(gradtemp,elem,hybrid,nets,nete) @@ -63,7 +68,7 @@ subroutine model_init2(elem,hybrid,deriv,hvcoord,tl,nets,nete ) else elem(ie)%state%w_i(:,:,nlevp,tl%n0) = (& elem(ie)%state%v(:,:,1,nlev,tl%n0)*elem(ie)%derived%gradphis(:,:,1) + & - elem(ie)%state%v(:,:,2,nlev,tl%n0)*elem(ie)%derived%gradphis(:,:,2))/g + elem(ie)%state%v(:,:,2,nlev,tl%n0)*elem(ie)%derived%gradphis(:,:,2))/g_from_phi(elem(ie)%state%phis(:,:)) ! DA_CHANGE endif ! assign phinh_i(nlevp) to be phis at all timelevels diff --git a/components/homme/src/theta-l/share/prim_advance_mod.F90 b/components/homme/src/theta-l/share/prim_advance_mod.F90 index fed34296c2bc..480e1619ad95 100644 --- a/components/homme/src/theta-l/share/prim_advance_mod.F90 +++ b/components/homme/src/theta-l/share/prim_advance_mod.F90 @@ -5,7 +5,7 @@ ! ! Man dynamics routines for "theta" nonhydrostatic model ! Original version: Mark Taylor 2017/1 -! +! call ! 2018/8 TOM sponge layer scaling from P. Lauritzen ! 09/2018: O. Guba code for new ftypes ! 2018/12: M. Taylor apply forcing assuming nearly constant p @@ -34,7 +34,7 @@ module prim_advance_mod use kinds, only: iulog, real_kind use perf_mod, only: t_adj_detailf, t_barrierf, t_startf, t_stopf ! _EXTERNAL use parallel_mod, only: abortmp, global_shared_buf, global_shared_sum, iam, parallel_t - use physical_constants, only: Cp, cp, cpwater_vapor, g, kappa, Rgas, Rwater_vapor, p0, TREF + use physical_constants, only: Cp, cp, cpwater_vapor, gravit, kappa, Rgas, Rwater_vapor, p0, TREF, rearth, omega_earth => omega use physics_mod, only: virtual_specific_heat, virtual_temperature use prim_si_mod, only: preq_vertadv_v1 use reduction_mod, only: parallelmax, reductionbuffer_ordered_1d_t @@ -82,6 +82,7 @@ end subroutine prim_advance_init1 !_____________________________________________________________________ subroutine prim_advance_exp(elem, deriv, hvcoord, hybrid,dt, tl, nets, nete, compute_diagnostics) use imex_mod, only: compute_stage_value_dirk + use deep_atm_mod, only: g_from_phi #ifdef ARKODE use arkode_mod, only: parameter_list, evolve_solution, & calc_nonlinear_stats, update_nonlinear_stats, & @@ -104,6 +105,7 @@ subroutine prim_advance_exp(elem, deriv, hvcoord, hybrid,dt, tl, nets, nete, co real (kind=real_kind) :: itertol,a1,a2,a3,a4,a5,a6,ahat1,ahat2 real (kind=real_kind) :: ahat3,ahat4,ahat5,ahat6,dhat1,dhat2,dhat3,dhat4 real (kind=real_kind) :: gamma,delta,ap,aphat,dhat5,offcenter + real (kind=real_kind) :: pnh(np, np, nlev), exner(np, np, nlev), dpnh_i(np,np,nlevp) integer :: ie,nm1,n0,np1,nstep,qsplit_stage,k integer :: n,i,j,maxiter @@ -141,7 +143,7 @@ subroutine prim_advance_exp(elem, deriv, hvcoord, hybrid,dt, tl, nets, nete, co ! this should not be needed, but in case physics update u without updating w b.c.: do ie=nets,nete elem(ie)%state%w_i(:,:,nlevp,n0) = (elem(ie)%state%v(:,:,1,nlev,n0)*elem(ie)%derived%gradphis(:,:,1) + & - elem(ie)%state%v(:,:,2,nlev,n0)*elem(ie)%derived%gradphis(:,:,2))/g + elem(ie)%state%v(:,:,2,nlev,n0)*elem(ie)%derived%gradphis(:,:,2))/g_from_phi(elem(ie)%state%phis(:,:)) !DA_CHANGE enddo #if !defined(CAM) && !defined(SCREAM) @@ -153,6 +155,7 @@ subroutine prim_advance_exp(elem, deriv, hvcoord, hybrid,dt, tl, nets, nete, co endif #endif + ! ================================== ! Take timestep ! ================================== @@ -517,7 +520,7 @@ end subroutine prim_advance_exp !----------------------------- APPLYCAMFORCING-DYNAMICS ---------------------------- subroutine applyCAMforcing_dynamics(elem,hvcoord,np1,dt,nets,nete) - + use deep_atm_mod, only: g_from_phi type (element_t) , intent(inout) :: elem(:) real (kind=real_kind), intent(in) :: dt type (hvcoord_t), intent(in) :: hvcoord @@ -537,7 +540,7 @@ subroutine applyCAMforcing_dynamics(elem,hvcoord,np1,dt,nets,nete) ! finally update w at the surface: elem(ie)%state%w_i(:,:,nlevp,np1) = (elem(ie)%state%v(:,:,1,nlev,np1)*elem(ie)%derived%gradphis(:,:,1) + & - elem(ie)%state%v(:,:,2,nlev,np1)*elem(ie)%derived%gradphis(:,:,2))/g + elem(ie)%state%v(:,:,2,nlev,np1)*elem(ie)%derived%gradphis(:,:,2))/g_from_phi(elem(ie)%state%phis(:,:)) !DA_CHANGE enddo end subroutine applyCAMforcing_dynamics @@ -546,6 +549,7 @@ end subroutine applyCAMforcing_dynamics !----------------------------- ADVANCE-HYPERVIS ---------------------------- subroutine advance_hypervis(elem,hvcoord,hybrid,deriv,nt,nets,nete,dt2,eta_ave_w) + use deep_atm_mod, only: g_from_phi, quasi_hydrostatic_terms ! ! take one timestep of: ! u(:,:,:,np) = u(:,:,:,np) + dt2*nu*laplacian**order ( u ) @@ -770,7 +774,7 @@ subroutine advance_hypervis(elem,hvcoord,hybrid,deriv,nt,nets,nete,dt2,eta_ave_w ! finally update w at the surface: elem(ie)%state%w_i(:,:,nlevp,nt) = (elem(ie)%state%v(:,:,1,nlev,nt)*elem(ie)%derived%gradphis(:,:,1) + & - elem(ie)%state%v(:,:,2,nlev,nt)*elem(ie)%derived%gradphis(:,:,2))/g + elem(ie)%state%v(:,:,2,nlev,nt)*elem(ie)%derived%gradphis(:,:,2))/g_from_phi(elem(ie)%state%phis(:,:)) ! DA_CHANGE enddo @@ -866,6 +870,7 @@ end subroutine advance_hypervis subroutine advance_physical_vis(elem,hvcoord,hybrid,deriv,nt,nets,nete,dt,mu_s,mu) + use deep_atm_mod, only: g_from_phi ! ! take one timestep of of physical viscosity (single laplace operator) for ! all state variables in both horizontal and vertical @@ -1006,7 +1011,7 @@ subroutine advance_physical_vis(elem,hvcoord,hybrid,deriv,nt,nets,nete,dt,mu_s,m ! finally update w at the surface: elem(ie)%state%w_i(:,:,nlevp,nt) = (elem(ie)%state%v(:,:,1,nlev,nt)*elem(ie)%derived%gradphis(:,:,1) + & - elem(ie)%state%v(:,:,2,nlev,nt)*elem(ie)%derived%gradphis(:,:,2))/g + elem(ie)%state%v(:,:,2,nlev,nt)*elem(ie)%derived%gradphis(:,:,2))/g_from_phi(elem(ie)%state%phis(:,:)) !DA_CHANGE enddo @@ -1022,6 +1027,9 @@ end subroutine advance_physical_vis subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,& deriv,nets,nete,compute_diagnostics,eta_ave_w,scale1,scale2,scale3) + use deep_atm_mod, only: r_hat_from_phi, z_from_phi, g_from_phi, quasi_hydrostatic_terms + use baroclinic_wave, only: baroclinic_wave_test + use control_mod, only: atm_is_deep ! =================================== ! compute the RHS, accumulate into u(np1) and apply DSS ! @@ -1088,7 +1096,18 @@ subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,& real (kind=real_kind) :: w_tens(np,np,nlevp) ! need to update w at surface as well real (kind=real_kind) :: theta_tens(np,np,nlev) real (kind=real_kind) :: phi_tens(np,np,nlevp) - + + real (kind=real_kind) :: recip_r_hat(np,np,nlev) + real (kind=real_kind) :: recip_r_hat_vec(np,np,2,nlev) + real (kind=real_kind) :: recip_r_hat_i_vec(np,np,2,nlevp) + real (kind=real_kind) :: recip_r_hat_i(np,np,nlevp) + real (kind=real_kind) :: fcos(np, np) + real (kind=real_kind) :: w_deep_metric_term(np,np,nlevp) + real (kind=real_kind) :: w_deep_coriolis_term(np,np,nlevp) + + real (kind=real_kind) :: u_deep_metric_term(np,np,nlev) + real (kind=real_kind) :: v_deep_metric_term(np,np,nlev) + real (kind=real_kind) :: u_deep_coriolis_term(np,np,nlev) real (kind=real_kind) :: pi(np,np,nlev) ! hydrostatic pressure real (kind=real_kind) :: pi_i(np,np,nlevp) ! hydrostatic pressure interfaces @@ -1098,6 +1117,7 @@ subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,& real (kind=real_kind) :: vtemp(np,np,2,nlev) ! generic gradient storage real (kind=real_kind), dimension(np,np) :: sdot_sum ! temporary field real (kind=real_kind) :: v1,v2,w,d_eta_dot_dpdn_dn, T0 + real (kind=real_kind) :: u_tmp, v_tmp, t_tmp, p_tmp, phis_tmp, rho_tmp, theta_tmp, q_tmp, psv_tmp, lon_tmp, lat_tmp, z_tmp, maxdiff integer :: i,j,k,kptr,ie, nlyr_tot call t_startf('compute_andor_apply_rhs') @@ -1114,43 +1134,42 @@ subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,& vtheta(:,:,:) = vtheta_dp(:,:,:)/dp3d(:,:,:) phi_i => elem(ie)%state%phinh_i(:,:,:,n0) -#ifdef ENERGY_DIAGNOSTICS - if (.not. theta_hydrostatic_mode) then - ! check w b.c. - temp(:,:,1) = (elem(ie)%state%v(:,:,1,nlev,n0)*elem(ie)%derived%gradphis(:,:,1) + & - elem(ie)%state%v(:,:,2,nlev,n0)*elem(ie)%derived%gradphis(:,:,2))/g - do j=1,np - do i=1,np - if ( abs(temp(i,j,1)-elem(ie)%state%w_i(i,j,nlevp,n0)) >1e-10) then - write(iulog,*) 'WARNING:CAAR w(n0) does not satisfy b.c.',ie,i,j,k - write(iulog,*) 'val1 = ',temp(i,j,1) - write(iulog,*) 'val2 = ',elem(ie)%state%w_i(i,j,nlevp,n0) - write(iulog,*) 'diff: ',temp(i,j,1)-elem(ie)%state%w_i(i,j,nlevp,n0) - endif - enddo - enddo - ! w boundary condition. just in case: - !elem(ie)%state%w_i(:,:,nlevp,n0) = (elem(ie)%state%v(:,:,1,nlev,n0)*elem(ie)%derived%gradphis(:,:,1) + & - ! elem(ie)%state%v(:,:,2,nlev,n0)*elem(ie)%derived%gradphis(:,:,2))/g - - ! check for layer spacing <= 1m - do k=1,nlev - do j=1,np - do i=1,np - if ((phi_i(i,j,k)-phi_i(i,j,k+1)) < g) then - write(iulog,*) 'WARNING:CAAR before ADV, delta z < 1m. ie,i,j,k=',ie,i,j,k - write(iulog,*) 'phi(i,j,k)= ',phi_i(i,j,k) - write(iulog,*) 'phi(i,j,k+1)=',phi_i(i,j,k+1) - endif - enddo - enddo - enddo - endif -#endif +!#ifdef ENERGY_DIAGNOSTICS +! if (.not. theta_hydrostatic_mode) then +! ! check w b.c. +! temp(:,:,1) = (elem(ie)%state%v(:,:,1,nlev,n0)*elem(ie)%derived%gradphis(:,:,1) + & +! elem(ie)%state%v(:,:,2,nlev,n0)*elem(ie)%derived%gradphis(:,:,2))/g +! do j=1,np +! do i=1,np +! if ( abs(temp(i,j,1)-elem(ie)%state%w_i(i,j,nlevp,n0)) >1e-10) then +! write(iulog,*) 'WARNING:CAAR w(n0) does not satisfy b.c.',ie,i,j,k +! write(iulog,*) 'val1 = ',temp(i,j,1) +! write(iulog,*) 'val2 = ',elem(ie)%state%w_i(i,j,nlevp,n0) +! write(iulog,*) 'diff: ',temp(i,j,1)-elem(ie)%state%w_i(i,j,nlevp,n0) +! endif +! enddo +! enddo +! ! w boundary condition. just in case: +! !elem(ie)%state%w_i(:,:,nlevp,n0) = (elem(ie)%state%v(:,:,1,nlev,n0)*elem(ie)%derived%gradphis(:,:,1) + & +! ! elem(ie)%state%v(:,:,2,nlev,n0)*elem(ie)%derived%gradphis(:,:,2))/g +! +! ! check for layer spacing <= 1m +! do k=1,nlev +! do j=1,np +! do i=1,np +! if ((phi_i(i,j,k)-phi_i(i,j,k+1)) < g) then +! write(iulog,*) 'WARNING:CAAR before ADV, delta z < 1m. ie,i,j,k=',ie,i,j,k +! write(iulog,*) 'phi(i,j,k)= ',phi_i(i,j,k) +! write(iulog,*) 'phi(i,j,k+1)=',phi_i(i,j,k+1) +! endif +! enddo +! enddo +! enddo +! endif +!#endif ! this routine will set dpnh_dp_i(nlevp)=1 - a very good approximation, that will ! then be corrected below, after the DSS. call pnh_and_exner_from_eos(hvcoord,vtheta_dp,dp3d,phi_i,pnh,exner,dpnh_dp_i,caller='CAAR') - dp3d_i(:,:,1) = dp3d(:,:,1) dp3d_i(:,:,nlevp) = dp3d(:,:,nlev) do k=2,nlev @@ -1176,10 +1195,15 @@ subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,& ! set to zero so H and NH can share code and reduce if statements elem(ie)%state%w_i(:,:,:,n0)=0 endif + recip_r_hat_i = 1.0_real_kind/r_hat_from_phi(phi_i, nlevp) ! DA_CHANGE + recip_r_hat_i_vec(:,:,1,:) = recip_r_hat_i ! DA_CHANGE + recip_r_hat_i_vec(:,:,1,:) = recip_r_hat_i ! DA_CHANGE do k=1,nlev phi(:,:,k) = (phi_i(:,:,k)+phi_i(:,:,k+1))/2 ! for diagnostics - + recip_r_hat(:,:,k) = 1.0_real_kind/r_hat_from_phi(phi(:,:,k)) !DA_CHANGE + recip_r_hat_vec(:,:,1,k) = recip_r_hat(:,:,k) + recip_r_hat_vec(:,:,2,k) = recip_r_hat(:,:,k) ! ================================ ! Accumulate mean Vel_rho flux in vn0 ! ================================ @@ -1187,8 +1211,8 @@ subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,& vtemp(:,:,2,k) = elem(ie)%state%v(:,:,2,k,n0)*dp3d(:,:,k) elem(ie)%derived%vn0(:,:,:,k)=elem(ie)%derived%vn0(:,:,:,k)+eta_ave_w*vtemp(:,:,:,k) - divdp(:,:,k)=divergence_sphere(vtemp(:,:,:,k),deriv,elem(ie)) - vort(:,:,k)=vorticity_sphere(elem(ie)%state%v(:,:,:,k,n0),deriv,elem(ie)) + divdp(:,:,k)=divergence_sphere(recip_r_hat_vec(:,:,:,k)*vtemp(:,:,:,k),deriv,elem(ie)) + vort(:,:,k)=recip_r_hat(:,:,k)*vorticity_sphere(elem(ie)%state%v(:,:,:,k,n0),deriv,elem(ie)) enddo ! Compute omega = Dpi/Dt Used only as a DIAGNOSTIC @@ -1204,7 +1228,7 @@ subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,& #else pi(:,:,k)=pi_i(:,:,k) + dp3d(:,:,k)/2 #endif - vtemp(:,:,:,k) = gradient_sphere( pi(:,:,k), deriv, elem(ie)%Dinv); + vtemp(:,:,:,k) = recip_r_hat_vec(:,:,:,k)* gradient_sphere( pi(:,:,k), deriv, elem(ie)%Dinv); vgrad_p(:,:,k) = elem(ie)%state%v(:,:,1,k,n0)*vtemp(:,:,1,k)+& elem(ie)%state%v(:,:,2,k,n0)*vtemp(:,:,2,k) omega(:,:,k) = (vgrad_p(:,:,k) - ( omega_i(:,:,k)+omega_i(:,:,k+1))/2) @@ -1222,6 +1246,9 @@ subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,& theta_vadv=0 v_vadv=0 else + if (atm_is_deep) then + call abortmp('ERROR: deep atmosphere not implemented for eulerian configurations') + end if sdot_sum=0 do k=1,nlev ! ================================================== @@ -1312,26 +1339,40 @@ subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,& enddo elem(ie)%derived%eta_dot_dpdn(:,:,nlev+1) = & elem(ie)%derived%eta_dot_dpdn(:,:,nlev+1) + eta_ave_w*eta_dot_dpdn(:,:,nlev+1) - + do j=1,np + do i=1,np + fcos(i,j) = 2.0D0*omega_earth*cos(elem(ie)%spherep(i,j)%lat) + end do + end do + ! ================================================ ! w,phi tendencies including surface ! ================================================ do k=1,nlev ! compute gradphi at interfaces and then average to levels - gradphinh_i(:,:,:,k) = gradient_sphere(phi_i(:,:,k),deriv,elem(ie)%Dinv) - - gradw_i(:,:,:,k) = gradient_sphere(elem(ie)%state%w_i(:,:,k,n0),deriv,elem(ie)%Dinv) + gradphinh_i(:,:,:,k) = recip_r_hat_vec(:,:,:,k)*gradient_sphere(phi_i(:,:,k),deriv,elem(ie)%Dinv) + gradw_i(:,:,:,k) = recip_r_hat_vec(:,:,:,k)*gradient_sphere(elem(ie)%state%w_i(:,:,k,n0),deriv,elem(ie)%Dinv) v_gradw_i(:,:,k) = v_i(:,:,1,k)*gradw_i(:,:,1,k) + v_i(:,:,2,k)*gradw_i(:,:,2,k) + if ( atm_is_deep .and. (k /= 1)) then + w_deep_metric_term(:,:,k) = (v_i(:,:,1,k)**2 + v_i(:,:,2,k)**2)/(rearth + z_from_phi(phi_i(:,:,k))) ! DA_CHANGE + w_deep_coriolis_term(:,:,k) = fcos * v_i(:,:,1,k) + else + w_deep_metric_term(:, :, k) = 0.0_real_kind + w_deep_coriolis_term(:,:,k) = 0.0_real_kind + end if ! w - tendency on interfaces - w_tens(:,:,k) = (-w_vadv_i(:,:,k) - v_gradw_i(:,:,k))*scale1 - scale2*g*(1-dpnh_dp_i(:,:,k) ) + w_tens(:,:,k) = (-w_vadv_i(:,:,k) - v_gradw_i(:,:,k))*scale1 & + + scale1*(w_deep_metric_term(:,:,k) + w_deep_coriolis_term(:,:,k)) & + - scale2*g_from_phi((phi_i(:,:,k) + phi_i(:,:,k+1))/2.0)*(1-dpnh_dp_i(:,:,k) ) ! DA_CHANGE + !write(*,*) "max 1-mu at lev", k, ": ", maxval(abs(1- dpnh_dp_i(:,:,k))) ! phi - tendency on interfaces ! vtemp(:,:,:,k) = gradphinh_i(:,:,:,k) + & ! (scale2-1)*hvcoord%hybi(k)*elem(ie)%derived%gradphis(:,:,:) v_gradphinh_i(:,:,k) = v_i(:,:,1,k)*gradphinh_i(:,:,1,k) & +v_i(:,:,2,k)*gradphinh_i(:,:,2,k) phi_tens(:,:,k) = (-phi_vadv_i(:,:,k) - v_gradphinh_i(:,:,k))*scale1 & - + scale2*g*elem(ie)%state%w_i(:,:,k,n0) + + scale2*(g_from_phi((phi_i(:,:,k) + phi_i(:,:,k+1))/2.0))*elem(ie)%state%w_i(:,:,k,n0) ! DA_CHANGE if (scale1/=scale2) then ! add imex phi_h splitting ! use approximate phi_h = hybi*phis @@ -1345,18 +1386,22 @@ subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,& ! k =nlevp case, all terms in the imex methods are treated explicitly at the boundary k =nlevp + w_deep_metric_term(:,:,k) = (v_i(:,:,1,k)**2 + v_i(:,:,2,k)**2)/(rearth + z_from_phi(phi_i(:,:,k))) ! DA_CHANGE + + w_deep_coriolis_term(:,:,k) = fcos * v_i(:,:,1,k) ! compute gradphi at interfaces and then average to levels - gradphinh_i(:,:,:,k) = gradient_sphere(phi_i(:,:,k),deriv,elem(ie)%Dinv) - gradw_i(:,:,:,k) = gradient_sphere(elem(ie)%state%w_i(:,:,k,n0),deriv,elem(ie)%Dinv) + gradphinh_i(:,:,:,k) = recip_r_hat_vec(:,:,:,k)*gradient_sphere(phi_i(:,:,k),deriv,elem(ie)%Dinv) + gradw_i(:,:,:,k) = recip_r_hat_vec(:,:,:,k)*gradient_sphere(elem(ie)%state%w_i(:,:,k,n0),deriv,elem(ie)%Dinv) v_gradw_i(:,:,k) = v_i(:,:,1,k)*gradw_i(:,:,1,k) + v_i(:,:,2,k)*gradw_i(:,:,2,k) ! w - tendency on interfaces - w_tens(:,:,k) = (-w_vadv_i(:,:,k) - v_gradw_i(:,:,k))*scale1 - scale1*g*(1-dpnh_dp_i(:,:,k) ) + w_tens(:,:,k) = (-w_vadv_i(:,:,k) - v_gradw_i(:,:,k))*scale1 - scale1*(g_from_phi((phi_i(:,:,k) + phi_i(:,:,k+1)))/2.0)*(1-dpnh_dp_i(:,:,k) ) & !DA_CHANGE + - scale1*(w_deep_metric_term(:,:,k) + w_deep_coriolis_term(:,:,k)) ! phi - tendency on interfaces v_gradphinh_i(:,:,k) = v_i(:,:,1,k)*gradphinh_i(:,:,1,k) & +v_i(:,:,2,k)*gradphinh_i(:,:,2,k) phi_tens(:,:,k) = (-phi_vadv_i(:,:,k) - v_gradphinh_i(:,:,k))*scale1 & - + scale1*g*elem(ie)%state%w_i(:,:,k,n0) + + scale1*(g_from_phi((phi_i(:,:,k) + phi_i(:,:,k+1))/2.0))*elem(ie)%state%w_i(:,:,k,n0) !DA_CHANGE @@ -1370,10 +1415,10 @@ subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,& if (theta_advect_form==0) then v_theta(:,:,1,k)=elem(ie)%state%v(:,:,1,k,n0)*vtheta_dp(:,:,k) v_theta(:,:,2,k)=elem(ie)%state%v(:,:,2,k,n0)*vtheta_dp(:,:,k) - div_v_theta(:,:,k)=divergence_sphere(v_theta(:,:,:,k),deriv,elem(ie)) + div_v_theta(:,:,k)=divergence_sphere(recip_r_hat_vec(:,:,:,k) * v_theta(:,:,:,k),deriv,elem(ie)) else ! alternate form, non-conservative, better HS topography results - v_theta(:,:,:,k) = gradient_sphere(vtheta(:,:,k),deriv,elem(ie)%Dinv) + v_theta(:,:,:,k) = recip_r_hat_vec(:, :, :, k) * gradient_sphere(vtheta(:,:,k),deriv,elem(ie)%Dinv) div_v_theta(:,:,k)=vtheta(:,:,k)*divdp(:,:,k) + & dp3d(:,:,k)*elem(ie)%state%v(:,:,1,k,n0)*v_theta(:,:,1,k) + & dp3d(:,:,k)*elem(ie)%state%v(:,:,2,k,n0)*v_theta(:,:,2,k) @@ -1387,21 +1432,21 @@ subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,& ! w vorticity correction term temp(:,:,k) = (elem(ie)%state%w_i(:,:,k,n0)**2 + & elem(ie)%state%w_i(:,:,k+1,n0)**2)/4 - wvor(:,:,:,k) = gradient_sphere(temp(:,:,k),deriv,elem(ie)%Dinv) + wvor(:,:,:,k) = recip_r_hat_vec(:,:,:,k) * gradient_sphere(temp(:,:,k),deriv,elem(ie)%Dinv) wvor(:,:,1,k) = wvor(:,:,1,k) - (elem(ie)%state%w_i(:,:,k,n0)*gradw_i(:,:,1,k) +& elem(ie)%state%w_i(:,:,k+1,n0)*gradw_i(:,:,1,k+1))/2 wvor(:,:,2,k) = wvor(:,:,2,k) - (elem(ie)%state%w_i(:,:,k,n0)*gradw_i(:,:,2,k) +& elem(ie)%state%w_i(:,:,k+1,n0)*gradw_i(:,:,2,k+1))/2 KE(:,:,k) = ( elem(ie)%state%v(:,:,1,k,n0)**2 + elem(ie)%state%v(:,:,2,k,n0)**2)/2 - gradKE(:,:,:,k) = gradient_sphere(KE(:,:,k),deriv,elem(ie)%Dinv) - gradexner(:,:,:,k) = gradient_sphere(exner(:,:,k),deriv,elem(ie)%Dinv) + gradKE(:,:,:,k) = recip_r_hat_vec(:,:,:,k)*gradient_sphere(KE(:,:,k),deriv,elem(ie)%Dinv) + gradexner(:,:,:,k) = recip_r_hat_vec(:,:,:,k)*gradient_sphere(exner(:,:,k),deriv,elem(ie)%Dinv) #if 0 ! another form: (good results in dcmip2012 test2.0) max=0.195 ! but bad results with HS topo ! grad(exner) =( grad(theta*exner) - exner*grad(theta))/theta - vtemp(:,:,:,k) = gradient_sphere(vtheta(:,:,k)*exner(:,:,k),deriv,elem(ie)%Dinv) - v_theta(:,:,:,k) = gradient_sphere(vtheta(:,:,k),deriv,elem(ie)%Dinv) + vtemp(:,:,:,k) = recip_r_hat(:,:,k)*gradient_sphere(vtheta(:,:,k)*exner(:,:,k),deriv,elem(ie)%Dinv) + v_theta(:,:,:,k) = recip_r_hat(:,:,k)*gradient_sphere(vtheta(:,:,k),deriv,elem(ie)%Dinv) gradexner(:,:,1,k) = (vtemp(:,:,1,k)-exner(:,:,k)*v_theta(:,:,1,k))/& vtheta(:,:,k) gradexner(:,:,2,k) = (vtemp(:,:,2,k)-exner(:,:,k)*v_theta(:,:,2,k))/& @@ -1409,8 +1454,8 @@ subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,& #endif #if 0 ! entropy form: dcmip2012 test2.0 best: max=0.130 (0.124 with conservation form theta) - vtemp(:,:,:,k) = gradient_sphere(vtheta(:,:,k)*exner(:,:,k),deriv,elem(ie)%Dinv) - v_theta(:,:,:,k) = gradient_sphere(log(vtheta(:,:,k)),deriv,elem(ie)%Dinv) + vtemp(:,:,:,k) = recip_r_hat(:,:,k)*gradient_sphere(vtheta(:,:,k)*exner(:,:,k),deriv,elem(ie)%Dinv) + v_theta(:,:,:,k) = recip_r_hat(:,:,k)gradient_sphere(log(vtheta(:,:,k)),deriv,elem(ie)%Dinv) gradexner(:,:,1,k) = (vtemp(:,:,1,k)-exner(:,:,k)*vtheta(:,:,k)*v_theta(:,:,1,k))/& vtheta(:,:,k) gradexner(:,:,2,k) = (vtemp(:,:,2,k)-exner(:,:,k)*vtheta(:,:,k)*v_theta(:,:,2,k))/& @@ -1419,7 +1464,7 @@ subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,& #if 0 ! another form: terrible results in dcmip2012 test2.0 ! grad(exner) = grad(p) * kappa * exner / p - gradexner(:,:,:,k) = gradient_sphere(pnh(:,:,k),deriv,elem(ie)%Dinv) + gradexner(:,:,:,k) = recip_r_hat(:,:,k)*gradient_sphere(pnh(:,:,k),deriv,elem(ie)%Dinv) gradexner(:,:,1,k) = gradexner(:,:,1,k)*(Rgas/Cp)*exner(:,:,k)/pnh(:,:,k) gradexner(:,:,2,k) = gradexner(:,:,2,k)*(Rgas/Cp)*exner(:,:,k)/pnh(:,:,k) #endif @@ -1431,7 +1476,7 @@ subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,& dpnh_dp_i(:,:,k+1)*gradphinh_i(:,:,2,k+1))/2 if (pgrad_correction==1) then - T0 = TREF-tref_lapse_rate*TREF*Cp/g ! = 97 + T0 = TREF-tref_lapse_rate*TREF*Cp/gravit ! = 97 DA_CHANGE #ifdef HOMMEXX_BFB_TESTING ! For BFB testing, calculate log(exner) using cxx_log() ! and then call gradient sphere. @@ -1441,14 +1486,28 @@ subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,& end do end do - vtemp(:,:,:,k)=gradient_sphere(temp(:,:,k),deriv,elem(ie)%Dinv) + vtemp(:,:,:,k)=recip_r_hat(:,:,k)*gradient_sphere(temp(:,:,k),deriv,elem(ie)%Dinv) #else - vtemp(:,:,:,k)=gradient_sphere(log(exner(:,:,k)),deriv,elem(ie)%Dinv) + vtemp(:,:,:,k)=recip_r_hat_vec(:,:,:,k)*gradient_sphere(log(exner(:,:,k)),deriv,elem(ie)%Dinv) #endif mgrad(:,:,1,k)=mgrad(:,:,1,k) + Cp*T0*(vtemp(:,:,1,k)-gradexner(:,:,1,k)/exner(:,:,k)) mgrad(:,:,2,k)=mgrad(:,:,2,k) + Cp*T0*(vtemp(:,:,2,k)-gradexner(:,:,2,k)/exner(:,:,k)) endif + if (atm_is_deep) then + u_deep_metric_term(:,:,k) = elem(ie)%state%v(:,:,1,k,n0) * & + (elem(ie)%state%w_i(:,:,k,n0)+elem(ie)%state%w_i(:,:,k+1,n0))/2 & + /(rearth + z_from_phi(phi(:,:,k))) !DA_CHANGE + v_deep_metric_term(:,:,k) = elem(ie)%state%v(:,:,2,k,n0) * & + (elem(ie)%state%w_i(:,:,k,n0)+elem(ie)%state%w_i(:,:,k+1,n0))/2 & + /(rearth + z_from_phi(phi(:,:,k))) ! DA_CHANGE + u_deep_coriolis_term(:,:,k) = fcos * & + (elem(ie)%state%w_i(:,:,k,n0)+elem(ie)%state%w_i(:,:,k+1,n0))/2 + else + u_deep_metric_term(:,:,k) = 0 + v_deep_metric_term(:,:,k) = 0 + u_deep_coriolis_term(:,:,k) = 0 + end if do j=1,np do i=1,np @@ -1470,6 +1529,8 @@ subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,& + v2*(elem(ie)%fcor(i,j) + vort(i,j,k)) & - gradKE(i,j,1,k) - mgrad(i,j,1,k) & -Cp*vtheta(i,j,k)*gradexner(i,j,1,k)& + -u_deep_metric_term(i,j,k) & + -u_deep_coriolis_term(i,j,k) & -wvor(i,j,1,k) )*scale1 @@ -1477,6 +1538,7 @@ subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,& - v1*(elem(ie)%fcor(i,j) + vort(i,j,k)) & - gradKE(i,j,2,k) - mgrad(i,j,2,k) & -Cp*vtheta(i,j,k)*gradexner(i,j,2,k) & + -v_deep_metric_term(i,j,k) & -wvor(i,j,2,k) )*scale1 #endif end do @@ -1484,145 +1546,144 @@ subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,& end do - -#ifdef ENERGY_DIAGNOSTICS - ! ========================================================= - ! diagnostics. not performance critical, dont thread - ! ========================================================= - if (compute_diagnostics) then - elem(ie)%accum%KEu_horiz1=0 - elem(ie)%accum%KEu_horiz2=0 - elem(ie)%accum%KEu_vert1=0 - elem(ie)%accum%KEu_vert2=0 - elem(ie)%accum%KEw_horiz1=0 - elem(ie)%accum%KEw_horiz2=0 - elem(ie)%accum%KEw_horiz3=0 - elem(ie)%accum%KEw_vert1=0 - elem(ie)%accum%KEw_vert2=0 - - elem(ie)%accum%PEhoriz1=0 - elem(ie)%accum%PEhoriz2=0 - elem(ie)%accum%IEvert1=0 - elem(ie)%accum%IEvert2=0 - elem(ie)%accum%PEvert1=0 - elem(ie)%accum%PEvert2=0 - elem(ie)%accum%T01=0 - elem(ie)%accum%T2=0 - elem(ie)%accum%S1=0 - elem(ie)%accum%S2=0 - elem(ie)%accum%P1=0 - elem(ie)%accum%P2=0 - - do k =1,nlev - do j=1,np - do i=1,np - d_eta_dot_dpdn_dn=(eta_dot_dpdn(i,j,k+1)-eta_dot_dpdn(i,j,k)) - ! Form horiz advection of KE-u - elem(ie)%accum%KEu_horiz1(i,j)=elem(ie)%accum%KEu_horiz1(i,j) & - -dp3d(i,j,k)*( & - elem(ie)%state%v(i,j,1,k,n0)*gradKE(i,j,1,k) + & - elem(ie)%state%v(i,j,2,k,n0)*gradKE(i,j,2,k) ) - elem(ie)%accum%KEu_horiz2(i,j)=elem(ie)%accum%KEu_horiz2(i,j) & - -KE(i,j,k)*divdp(i,j,k) - ! Form horiz advection of KE-w - elem(ie)%accum%KEw_horiz1(i,j)=elem(ie)%accum%KEw_horiz1(i,j)- & - dp3d(i,j,k) * (& - elem(ie)%state%w_i(i,j,k,n0) * v_gradw_i(i,j,k) + & - elem(ie)%state%w_i(i,j,k+1,n0) * v_gradw_i(i,j,k+1) )/2 - elem(ie)%accum%KEw_horiz2(i,j)=elem(ie)%accum%KEw_horiz2(i,j)- & - divdp(i,j,k)*(elem(ie)%state%w_i(i,j,k,n0)**2 + & - elem(ie)%state%w_i(i,j,k+1,n0)**2 ) /4 - elem(ie)%accum%KEw_horiz3(i,j)=elem(ie)%accum%KEw_horiz3(i,j) & - -dp3d(i,j,k) * (elem(ie)%state%v(i,j,1,k,n0) * wvor(i,j,1,k) + & - elem(ie)%state%v(i,j,2,k,n0) * wvor(i,j,2,k)) - ! Form vertical advection of KE-u - elem(ie)%accum%KEu_vert1(i,j)=elem(ie)%accum%KEu_vert1(i,j)- & - (elem(ie)%state%v(i,j,1,k,n0) * v_vadv(i,j,1,k) + & - elem(ie)%state%v(i,j,2,k,n0) *v_vadv(i,j,2,k))*dp3d(i,j,k) - elem(ie)%accum%KEu_vert2(i,j)=elem(ie)%accum%KEu_vert2(i,j)- & - 0.5*((elem(ie)%state%v(i,j,1,k,n0))**2 + & - (elem(ie)%state%v(i,j,2,k,n0))**2)*d_eta_dot_dpdn_dn - ! Form vertical advection of KE-w - elem(ie)%accum%KEw_vert1(i,j)=elem(ie)%accum%KEw_vert1(i,j) - & - dp3d(i,j,k) * & - (w_vadv_i(i,j,k)*elem(ie)%state%w_i(i,j,k,n0)+ & - w_vadv_i(i,j,k+1)*elem(ie)%state%w_i(i,j,k+1,n0))/2 - - elem(ie)%accum%KEw_vert2(i,j)=elem(ie)%accum%KEw_vert2(i,j) & - -d_eta_dot_dpdn_dn* & - (.5*elem(ie)%state%w_i(i,j,k,n0)**2 +& - .5*elem(ie)%state%w_i(i,j,k+1,n0)**2)/2 - - ! Form IEvert1 - elem(ie)%accum%IEvert1(i,j)=elem(ie)%accum%IEvert1(i,j) & - -Cp*exner(i,j,k)*theta_vadv(i,j,k) - ! Form IEvert2 - ! here use of dpnh_dp_i on boundry (with incorrect data) - ! is harmess becuase eta_dot_dpdn=0 - elem(ie)%accum%IEvert2(i,j)=elem(ie)%accum%IEvert2(i,j) & - + ( dpnh_dp_i(i,j,k)*eta_dot_dpdn(i,j,k)+ & - dpnh_dp_i(i,j,k+1)*eta_dot_dpdn(i,j,k+1)) & - *(phi_i(i,j,k+1)-phi_i(i,j,k))/2 - - ! Form PEhoriz1 - elem(ie)%accum%PEhoriz1(i,j)=(elem(ie)%accum%PEhoriz1(i,j)) & - -phi(i,j,k)*divdp(i,j,k) - ! Form PEhoriz2 - elem(ie)%accum%PEhoriz2(i,j)=elem(ie)%accum%PEhoriz2(i,j) & - -dp3d(i,j,k)* & - (elem(ie)%state%v(i,j,1,k,n0)* & - (gradphinh_i(i,j,1,k)+gradphinh_i(i,j,1,k+1))/2 + & - elem(ie)%state%v(i,j,2,k,n0)* & - (gradphinh_i(i,j,2,k)+gradphinh_i(i,j,2,k+1))/2 ) - - ! Form PEvert1 - elem(ie)%accum%PEvert1(i,j) = elem(ie)%accum%PEvert1(i,j) & - -phi(i,j,k)*d_eta_dot_dpdn_dn - elem(ie)%accum%PEvert2(i,j) = elem(ie)%accum%PEvert2(i,j) & - -dp3d(i,j,k)*(phi_vadv_i(i,j,k)+phi_vadv_i(i,j,k+1))/2 - - ! Form T01 - elem(ie)%accum%T01(i,j)=elem(ie)%accum%T01(i,j) & - -(Cp*elem(ie)%state%vtheta_dp(i,j,k,n0)) & - *(gradexner(i,j,1,k)*elem(ie)%state%v(i,j,1,k,n0) + & - gradexner(i,j,2,k)*elem(ie)%state%v(i,j,2,k,n0)) - ! Form S1 - elem(ie)%accum%S1(i,j)=elem(ie)%accum%S1(i,j) & - -Cp*exner(i,j,k)*div_v_theta(i,j,k) - - ! Form P1 = -P2 (no reason to compute P2?) - elem(ie)%accum%P1(i,j)=elem(ie)%accum%P1(i,j) -g*dp3d(i,j,k)* & - ( elem(ie)%state%w_i(i,j,k,n0) + & - elem(ie)%state%w_i(i,j,k+1,n0) )/2 - ! Form P2 - elem(ie)%accum%P2(i,j)=elem(ie)%accum%P2(i,j) + g*dp3d(i,j,k)*& - ( elem(ie)%state%w_i(i,j,k,n0) + & - elem(ie)%state%w_i(i,j,k+1,n0) )/2 - enddo - enddo - enddo - - ! these terms are better easier to compute by summing interfaces - do k=2,nlev - elem(ie)%accum%T2(:,:)=elem(ie)%accum%T2(:,:)+ & - (g*elem(ie)%state%w_i(:,:,k,n0)-v_gradphinh_i(:,:,k)) & - * dpnh_dp_i(:,:,k)*dp3d_i(:,:,k) - enddo - ! boundary terms - do k=1,nlevp,nlev - elem(ie)%accum%T2(:,:)=elem(ie)%accum%T2(:,:)+ & - (g*elem(ie)%state%w_i(:,:,k,n0)-v_gradphinh_i(:,:,k)) & - * dpnh_dp_i(:,:,k)*dp3d_i(:,:,k)/2 - enddo - ! boundary term is incorrect. save the term so we can correct it - ! once we have coorect value of dpnh_dp_i: - elem(ie)%accum%T2_nlevp_term(:,:)=& - (g*elem(ie)%state%w_i(:,:,nlevp,n0)-v_gradphinh_i(:,:,nlevp)) & - * dp3d_i(:,:,nlevp)/2 - - endif -#endif +!#ifdef ENERGY_DIAGNOSTICS +! ! ========================================================= +! ! diagnostics. not performance critical, dont thread +! ! ========================================================= +! if (compute_diagnostics) then +! elem(ie)%accum%KEu_horiz1=0 +! elem(ie)%accum%KEu_horiz2=0 +! elem(ie)%accum%KEu_vert1=0 +! elem(ie)%accum%KEu_vert2=0 +! elem(ie)%accum%KEw_horiz1=0 +! elem(ie)%accum%KEw_horiz2=0 +! elem(ie)%accum%KEw_horiz3=0 +! elem(ie)%accum%KEw_vert1=0 +! elem(ie)%accum%KEw_vert2=0 +! +! elem(ie)%accum%PEhoriz1=0 +! elem(ie)%accum%PEhoriz2=0 +! elem(ie)%accum%IEvert1=0 +! elem(ie)%accum%IEvert2=0 +! elem(ie)%accum%PEvert1=0 +! elem(ie)%accum%PEvert2=0 +! elem(ie)%accum%T01=0 +! elem(ie)%accum%T2=0 +! elem(ie)%accum%S1=0 +! elem(ie)%accum%S2=0 +! elem(ie)%accum%P1=0 +! elem(ie)%accum%P2=0 +! +! do k =1,nlev +! do j=1,np +! do i=1,np +! d_eta_dot_dpdn_dn=(eta_dot_dpdn(i,j,k+1)-eta_dot_dpdn(i,j,k)) +! ! Form horiz advection of KE-u +! elem(ie)%accum%KEu_horiz1(i,j)=elem(ie)%accum%KEu_horiz1(i,j) & +! -dp3d(i,j,k)*( & +! elem(ie)%state%v(i,j,1,k,n0)*gradKE(i,j,1,k) + & +! elem(ie)%state%v(i,j,2,k,n0)*gradKE(i,j,2,k) ) +! elem(ie)%accum%KEu_horiz2(i,j)=elem(ie)%accum%KEu_horiz2(i,j) & +! -KE(i,j,k)*divdp(i,j,k) +! ! Form horiz advection of KE-w +! elem(ie)%accum%KEw_horiz1(i,j)=elem(ie)%accum%KEw_horiz1(i,j)- & +! dp3d(i,j,k) * (& +! elem(ie)%state%w_i(i,j,k,n0) * v_gradw_i(i,j,k) + & +! elem(ie)%state%w_i(i,j,k+1,n0) * v_gradw_i(i,j,k+1) )/2 +! elem(ie)%accum%KEw_horiz2(i,j)=elem(ie)%accum%KEw_horiz2(i,j)- & +! divdp(i,j,k)*(elem(ie)%state%w_i(i,j,k,n0)**2 + & +! elem(ie)%state%w_i(i,j,k+1,n0)**2 ) /4 +! elem(ie)%accum%KEw_horiz3(i,j)=elem(ie)%accum%KEw_horiz3(i,j) & +! -dp3d(i,j,k) * (elem(ie)%state%v(i,j,1,k,n0) * wvor(i,j,1,k) + & +! elem(ie)%state%v(i,j,2,k,n0) * wvor(i,j,2,k)) +! ! Form vertical advection of KE-u +! elem(ie)%accum%KEu_vert1(i,j)=elem(ie)%accum%KEu_vert1(i,j)- & +! (elem(ie)%state%v(i,j,1,k,n0) * v_vadv(i,j,1,k) + & +! elem(ie)%state%v(i,j,2,k,n0) *v_vadv(i,j,2,k))*dp3d(i,j,k) +! elem(ie)%accum%KEu_vert2(i,j)=elem(ie)%accum%KEu_vert2(i,j)- & +! 0.5*((elem(ie)%state%v(i,j,1,k,n0))**2 + & +! (elem(ie)%state%v(i,j,2,k,n0))**2)*d_eta_dot_dpdn_dn +! ! Form vertical advection of KE-w +! elem(ie)%accum%KEw_vert1(i,j)=elem(ie)%accum%KEw_vert1(i,j) - & +! dp3d(i,j,k) * & +! (w_vadv_i(i,j,k)*elem(ie)%state%w_i(i,j,k,n0)+ & +! w_vadv_i(i,j,k+1)*elem(ie)%state%w_i(i,j,k+1,n0))/2 +! +! elem(ie)%accum%KEw_vert2(i,j)=elem(ie)%accum%KEw_vert2(i,j) & +! -d_eta_dot_dpdn_dn* & +! (.5*elem(ie)%state%w_i(i,j,k,n0)**2 +& +! .5*elem(ie)%state%w_i(i,j,k+1,n0)**2)/2 +! +! ! Form IEvert1 +! elem(ie)%accum%IEvert1(i,j)=elem(ie)%accum%IEvert1(i,j) & +! -Cp*exner(i,j,k)*theta_vadv(i,j,k) +! ! Form IEvert2 +! ! here use of dpnh_dp_i on boundry (with incorrect data) +! ! is harmess becuase eta_dot_dpdn=0 +! elem(ie)%accum%IEvert2(i,j)=elem(ie)%accum%IEvert2(i,j) & +! + ( dpnh_dp_i(i,j,k)*eta_dot_dpdn(i,j,k)+ & +! dpnh_dp_i(i,j,k+1)*eta_dot_dpdn(i,j,k+1)) & +! *(phi_i(i,j,k+1)-phi_i(i,j,k))/2 +! +! ! Form PEhoriz1 +! elem(ie)%accum%PEhoriz1(i,j)=(elem(ie)%accum%PEhoriz1(i,j)) & +! -phi(i,j,k)*divdp(i,j,k) +! ! Form PEhoriz2 +! elem(ie)%accum%PEhoriz2(i,j)=elem(ie)%accum%PEhoriz2(i,j) & +! -dp3d(i,j,k)* & +! (elem(ie)%state%v(i,j,1,k,n0)* & +! (gradphinh_i(i,j,1,k)+gradphinh_i(i,j,1,k+1))/2 + & +! elem(ie)%state%v(i,j,2,k,n0)* & +! (gradphinh_i(i,j,2,k)+gradphinh_i(i,j,2,k+1))/2 ) +! +! ! Form PEvert1 +! elem(ie)%accum%PEvert1(i,j) = elem(ie)%accum%PEvert1(i,j) & +! -phi(i,j,k)*d_eta_dot_dpdn_dn +! elem(ie)%accum%PEvert2(i,j) = elem(ie)%accum%PEvert2(i,j) & +! -dp3d(i,j,k)*(phi_vadv_i(i,j,k)+phi_vadv_i(i,j,k+1))/2 +! +! ! Form T01 +! elem(ie)%accum%T01(i,j)=elem(ie)%accum%T01(i,j) & +! -(Cp*elem(ie)%state%vtheta_dp(i,j,k,n0)) & +! *(gradexner(i,j,1,k)*elem(ie)%state%v(i,j,1,k,n0) + & +! gradexner(i,j,2,k)*elem(ie)%state%v(i,j,2,k,n0)) +! ! Form S1 +! elem(ie)%accum%S1(i,j)=elem(ie)%accum%S1(i,j) & +! -Cp*exner(i,j,k)*div_v_theta(i,j,k) +! +! ! Form P1 = -P2 (no reason to compute P2?) +! elem(ie)%accum%P1(i,j)=elem(ie)%accum%P1(i,j) -gravit*dp3d(i,j,k)* & !DA_CHANGE +! ( elem(ie)%state%w_i(i,j,k,n0) + & +! elem(ie)%state%w_i(i,j,k+1,n0) )/2 +! ! Form P2 +! elem(ie)%accum%P2(i,j)=elem(ie)%accum%P2(i,j) + gravit*dp3d(i,j,k)*& !DA_CHANGE +! ( elem(ie)%state%w_i(i,j,k,n0) + & +! elem(ie)%state%w_i(i,j,k+1,n0) )/2 +! enddo +! enddo +! enddo +! +! ! these terms are better easier to compute by summing interfaces +! do k=2,nlev +! elem(ie)%accum%T2(:,:)=elem(ie)%accum%T2(:,:)+ & +! (gravit*elem(ie)%state%w_i(:,:,k,n0)-v_gradphinh_i(:,:,k)) & !DA_CHANGE +! * dpnh_dp_i(:,:,k)*dp3d_i(:,:,k) +! enddo +! ! boundary terms +! do k=1,nlevp,nlev +! elem(ie)%accum%T2(:,:)=elem(ie)%accum%T2(:,:)+ & +! (gravit*elem(ie)%state%w_i(:,:,k,n0)-v_gradphinh_i(:,:,k)) & ! DA_CHANGE +! * dpnh_dp_i(:,:,k)*dp3d_i(:,:,k)/2 +! enddo +! ! boundary term is incorrect. save the term so we can correct it +! ! once we have coorect value of dpnh_dp_i: +! elem(ie)%accum%T2_nlevp_term(:,:)=& +! (gravit*elem(ie)%state%w_i(:,:,nlevp,n0)-v_gradphinh_i(:,:,nlevp)) & !DA_CHANGE +! * dp3d_i(:,:,nlevp)/2 +! +! endif +!#endif do k=1,nlev @@ -1732,56 +1793,56 @@ subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,& ! solve for (dpnh_dp_i-1) dpnh_dp_i(:,:,nlevp) = 1 + ( & ((elem(ie)%state%v(:,:,1,nlev,np1)*elem(ie)%derived%gradphis(:,:,1) + & - elem(ie)%state%v(:,:,2,nlev,np1)*elem(ie)%derived%gradphis(:,:,2))/g - & + elem(ie)%state%v(:,:,2,nlev,np1)*elem(ie)%derived%gradphis(:,:,2))/g_from_phi(elem(ie)%state%phis(:,:)) - & ! DA_CHANGE elem(ie)%state%w_i(:,:,nlevp,np1)) / & - (g + ( elem(ie)%derived%gradphis(:,:,1)**2 + & - elem(ie)%derived%gradphis(:,:,2)**2)/(2*g)) ) / dt2 + (g_from_phi(elem(ie)%state%phis(:,:)) + ( elem(ie)%derived%gradphis(:,:,1)**2 + & + elem(ie)%derived%gradphis(:,:,2)**2)/(2*g_from_phi(elem(ie)%state%phis(:,:))) ) )/ dt2 !DA_CHANGE ! update solution with new dpnh_dp_i value: elem(ie)%state%w_i(:,:,nlevp,np1) = elem(ie)%state%w_i(:,:,nlevp,np1) +& - scale1*dt2*g*(dpnh_dp_i(:,:,nlevp)-1) + scale1*dt2*g_from_phi(elem(ie)%state%phis(:,:))*(dpnh_dp_i(:,:,nlevp)-1) !DA_CHANGE elem(ie)%state%v(:,:,1,nlev,np1) = elem(ie)%state%v(:,:,1,nlev,np1) -& scale1*dt2*(dpnh_dp_i(:,:,nlevp)-1)*elem(ie)%derived%gradphis(:,:,1)/2 elem(ie)%state%v(:,:,2,nlev,np1) = elem(ie)%state%v(:,:,2,nlev,np1) -& scale1*dt2*(dpnh_dp_i(:,:,nlevp)-1)*elem(ie)%derived%gradphis(:,:,2)/2 -#ifdef ENERGY_DIAGNOSTICS - ! add in boundary term to T2 and S2 diagnostics: - if (compute_diagnostics) then - elem(ie)%accum%T2(:,:)=elem(ie)%accum%T2(:,:)+ & - elem(ie)%accum%T2_nlevp_term(:,:)*(dpnh_dp_i(:,:,nlevp)-1) - elem(ie)%accum%S2(:,:)=-elem(ie)%accum%T2(:,:) - endif - - ! check w b.c. - temp(:,:,1) = (elem(ie)%state%v(:,:,1,nlev,np1)*elem(ie)%derived%gradphis(:,:,1) + & - elem(ie)%state%v(:,:,2,nlev,np1)*elem(ie)%derived%gradphis(:,:,2))/g - do j=1,np - do i=1,np - if ( abs(temp(i,j,1)-elem(ie)%state%w_i(i,j,nlevp,np1)) >1e-10) then - write(iulog,*) 'WARNING:CAAR w(np1) does not satisfy b.c.',ie,i,j,k - write(iulog,*) 'val1 = ',temp(i,j,1) - write(iulog,*) 'val2 = ',elem(ie)%state%w_i(i,j,nlevp,np1) - write(iulog,*) 'diff: ',temp(i,j,1)-elem(ie)%state%w_i(i,j,nlevp,np1) - endif - enddo - enddo - - ! check for layer spacing <= 1m - if (scale3 /= 0) then - do k=1,nlev - do j=1,np - do i=1,np - if ((elem(ie)%state%phinh_i(i,j,k,np1)-elem(ie)%state%phinh_i(i,j,k+1,np1)) < g) then - write(iulog,*) 'WARNING:CAAR after ADV, delta z < 1m. ie,i,j,k=',ie,i,j,k - write(iulog,*) 'phi(i,j,k)= ',elem(ie)%state%phinh_i(i,j,k,np1) - write(iulog,*) 'phi(i,j,k+1)=',elem(ie)%state%phinh_i(i,j,k+1,np1) - endif - enddo - enddo - enddo - endif -#endif +!#ifdef ENERGY_DIAGNOSTICS +! ! add in boundary term to T2 and S2 diagnostics: +! if (compute_diagnostics) then +! elem(ie)%accum%T2(:,:)=elem(ie)%accum%T2(:,:)+ & +! elem(ie)%accum%T2_nlevp_term(:,:)*(dpnh_dp_i(:,:,nlevp)-1) +! elem(ie)%accum%S2(:,:)=-elem(ie)%accum%T2(:,:) +! endif +! +! ! check w b.c. +! temp(:,:,1) = (elem(ie)%state%v(:,:,1,nlev,np1)*elem(ie)%derived%gradphis(:,:,1) + & +! elem(ie)%state%v(:,:,2,nlev,np1)*elem(ie)%derived%gradphis(:,:,2))/gravit ! DA_CHANGE +! do j=1,np +! do i=1,np +! if ( abs(temp(i,j,1)-elem(ie)%state%w_i(i,j,nlevp,np1)) >1e-10) then +! write(iulog,*) 'WARNING:CAAR w(np1) does not satisfy b.c.',ie,i,j,k +! write(iulog,*) 'val1 = ',temp(i,j,1) +! write(iulog,*) 'val2 = ',elem(ie)%state%w_i(i,j,nlevp,np1) +! write(iulog,*) 'diff: ',temp(i,j,1)-elem(ie)%state%w_i(i,j,nlevp,np1) +! endif +! enddo +! enddo +! +! ! check for layer spacing <= 1m +! if (scale3 /= 0) then +! do k=1,nlev +! do j=1,np +! do i=1,np +! if ((elem(ie)%state%phinh_i(i,j,k,np1)-elem(ie)%state%phinh_i(i,j,k+1,np1)) < gravit) then !DA_CHANGE +! write(iulog,*) 'WARNING:CAAR after ADV, delta z < 1m. ie,i,j,k=',ie,i,j,k +! write(iulog,*) 'phi(i,j,k)= ',elem(ie)%state%phinh_i(i,j,k,np1) +! write(iulog,*) 'phi(i,j,k+1)=',elem(ie)%state%phinh_i(i,j,k+1,np1) +! endif +! enddo +! enddo +! enddo +! endif +!#endif endif if (scale3 /= 0) then call limiter_dp3d_k(elem(ie)%state%dp3d(:,:,:,np1),elem(ie)%state%vtheta_dp(:,:,:,np1),& diff --git a/components/homme/src/theta-l/share/prim_state_mod.F90 b/components/homme/src/theta-l/share/prim_state_mod.F90 index 5a571701e584..3dcd8e93ad18 100644 --- a/components/homme/src/theta-l/share/prim_state_mod.F90 +++ b/components/homme/src/theta-l/share/prim_state_mod.F90 @@ -25,7 +25,7 @@ module prim_state_mod use viscosity_mod, only: compute_zeta_C0 use reduction_mod, only: parallelmax,parallelmin,parallelmaxwithindex,parallelminwithindex use perf_mod, only: t_startf, t_stopf - use physical_constants, only : p0,Cp,g,Rgas + use physical_constants, only : p0,Cp,gravit,Rgas implicit none private @@ -80,6 +80,7 @@ end subroutine prim_printstate_init subroutine prim_printstate(elem, tl,hybrid,hvcoord,nets,nete) use physical_constants, only: dd_pi + use deep_atm_mod, only: g_from_phi type(element_t), intent(inout), target :: elem(:) type(TimeLevel_t),target, intent(in) :: tl @@ -240,8 +241,7 @@ subroutine prim_printstate(elem, tl,hybrid,hvcoord,nets,nete) do k=1,nlev dphi(:,:,k)=-(phi_i(:,:,k+1)-phi_i(:,:,k)) w_over_dz(:,:,k)=& - g*max(elem(ie)%state%w_i(:,:,k,n0)/dphi(:,:,k),& - elem(ie)%state%w_i(:,:,k+1,n0)/dphi(:,:,k)) + g_from_phi(phi_i(:,:,k))*max(elem(ie)%state%w_i(:,:,k,n0)/dphi(:,:,k),elem(ie)%state%w_i(:,:,k+1,n0)/dphi(:,:,k)) ! DA_CHANGE enddo ! surface pressure @@ -368,7 +368,7 @@ subroutine prim_printstate(elem, tl,hybrid,hvcoord,nets,nete) thetasum_p= global_shared_sum(12) - scale=1/g ! assume code is using Pa + scale=1/gravit ! assume code is using Pa WARNING: not adapted to Deep Atmosphere if (hvcoord%ps0 < 2000 ) scale=100*scale ! code is using mb ! after scaling, Energy is in J/m**2, Mass kg/m**2 ! so rate of changes are W/m**2 @@ -404,8 +404,8 @@ subroutine prim_printstate(elem, tl,hybrid,hvcoord,nets,nete) write(iulog,109) "mu = ",mumin_local(1)," (",nint(mumin_local(2)),")",& mumax_local(1)," (",nint(mumax_local(2)),")" - write(iulog,109) "dz(m) = ",phimin_local(1)/g," (",nint(phimin_local(2)),")",& - phimax_local(1)/g," (",nint(phimax_local(2)),")",phisum_p/g + write(iulog,109) "dz(m) = ",phimin_local(1)/gravit," (",nint(phimin_local(2)),")",& + phimax_local(1)/gravit," (",nint(phimax_local(2)),")",phisum_p/gravit !DA_CHANGE write(iulog,109) "dp = ",dpmin_local(1)," (",nint(dpmin_local(2)),")",& dpmax_local(1)," (",nint(dpmax_local(2)),")",dpsum_p diff --git a/components/homme/src/theta-l/share/vertremap_mod.F90 b/components/homme/src/theta-l/share/vertremap_mod.F90 index b14a413d64bc..cd98b29d5794 100644 --- a/components/homme/src/theta-l/share/vertremap_mod.F90 +++ b/components/homme/src/theta-l/share/vertremap_mod.F90 @@ -38,7 +38,8 @@ subroutine vertical_remap(hybrid,elem,hvcoord,dt,np1,np1_qdp,nets,nete) use hybvcoord_mod, only: hvcoord_t use control_mod, only: rsplit use hybrid_mod, only: hybrid_t - use physical_constants, only : Cp,g + use physical_constants, only: gravit + use deep_atm_mod, only: g_from_phi, r_hat_from_phi type (hybrid_t), intent(in) :: hybrid ! distributed parallel structure (shared) type (element_t), intent(inout) :: elem(:) @@ -48,6 +49,7 @@ subroutine vertical_remap(hybrid,elem,hvcoord,dt,np1,np1_qdp,nets,nete) integer :: ie,i,j,k,np1,nets,nete,np1_qdp integer :: q + real (kind=real_kind), dimension(np,np) :: r_hat real (kind=real_kind), dimension(np,np,nlev) :: dp,dp_star real (kind=real_kind), dimension(np,np,nlevp) :: phi_ref real (kind=real_kind), dimension(np,np,nlev,5) :: ttmp @@ -78,8 +80,8 @@ subroutine vertical_remap(hybrid,elem,hvcoord,dt,np1,np1_qdp,nets,nete) elem(ie)%state%ps_v(:,:,np1) = hvcoord%hyai(1)*hvcoord%ps0 + & sum(elem(ie)%state%dp3d(:,:,:,np1),3) do k=1,nlev - dp(:,:,k) = ( hvcoord%hyai(k+1) - hvcoord%hyai(k) )*hvcoord%ps0 + & - ( hvcoord%hybi(k+1) - hvcoord%hybi(k) )*elem(ie)%state%ps_v(:,:,np1) + dp(:,:,k) = (( hvcoord%hyai(k+1) - hvcoord%hyai(k) )*hvcoord%ps0 + & + ( hvcoord%hybi(k+1) - hvcoord%hybi(k) )*elem(ie)%state%ps_v(:,:,np1)) if (rsplit==0) then dp_star(:,:,k) = dp(:,:,k) + dt*(elem(ie)%derived%eta_dot_dpdn(:,:,k+1) -& elem(ie)%derived%eta_dot_dpdn(:,:,k)) @@ -105,7 +107,7 @@ subroutine vertical_remap(hybrid,elem,hvcoord,dt,np1,np1_qdp,nets,nete) if (rsplit>0) then ! remove hydrostatic phi befor remap - call phi_from_eos(hvcoord,elem(ie)%state%phis,elem(ie)%state%vtheta_dp(:,:,:,np1),dp_star,phi_ref) + call phi_from_eos(hvcoord,elem(ie)%state%phis, elem(ie)%state%ps_v(:,:,np1),elem(ie)%state%vtheta_dp(:,:,:,np1),dp_star,phi_ref) elem(ie)%state%phinh_i(:,:,:,np1)=& elem(ie)%state%phinh_i(:,:,:,np1) -phi_ref(:,:,:) @@ -139,14 +141,14 @@ subroutine vertical_remap(hybrid,elem,hvcoord,dt,np1,np1_qdp,nets,nete) enddo ! depends on theta, so do this after updating theta: - call phi_from_eos(hvcoord,elem(ie)%state%phis,elem(ie)%state%vtheta_dp(:,:,:,np1),dp,phi_ref) + call phi_from_eos(hvcoord,elem(ie)%state%phis,elem(ie)%state%ps_v(:,:,np1),elem(ie)%state%vtheta_dp(:,:,:,np1),dp,phi_ref) elem(ie)%state%phinh_i(:,:,:,np1)=& elem(ie)%state%phinh_i(:,:,:,np1)+phi_ref(:,:,:) ! since u changed, update w b.c.: elem(ie)%state%w_i(:,:,nlevp,np1) = (elem(ie)%state%v(:,:,1,nlev,np1)*elem(ie)%derived%gradphis(:,:,1) + & - elem(ie)%state%v(:,:,2,nlev,np1)*elem(ie)%derived%gradphis(:,:,2))/g + elem(ie)%state%v(:,:,2,nlev,np1)*elem(ie)%derived%gradphis(:,:,2))/g_from_phi(elem(ie)%state%phis(:,:)) !DA_CHANGE endif