diff --git a/bld/build-namelist b/bld/build-namelist index 8408af37bb..ad43b5eabe 100755 --- a/bld/build-namelist +++ b/bld/build-namelist @@ -3619,6 +3619,9 @@ if ($dyn =~ /se/) { se_nu_div se_nu_p se_nu_top + se_sponge_del4_nu_fac + se_sponge_del4_nu_div_fac + se_sponge_del4_lev se_qsplit se_rsplit se_statediag_numtrac @@ -3634,9 +3637,6 @@ if ($dyn =~ /se/) { se_kmin_jet se_kmax_jet se_phys_dyn_cp - se_raytau0 - se_raykrange - se_rayk0 se_molecular_diff ); diff --git a/bld/namelist_files/namelist_defaults_cam.xml b/bld/namelist_files/namelist_defaults_cam.xml index f511cbcf1e..c855b80a82 100644 --- a/bld/namelist_files/namelist_defaults_cam.xml +++ b/bld/namelist_files/namelist_defaults_cam.xml @@ -649,11 +649,11 @@ 18 -0.1D0 +0.7D0 0.4D0 0.55D0 0.5D0 -0.325D0 +0.5D0 0.7D0 0.5D0 0.5D0 @@ -1996,6 +1996,11 @@ 1 1 3 + + 2 + 1 + 2 + 1 1.0D0 @@ -2658,9 +2663,6 @@ - 0 - 240 - 2 .false. @@ -2668,60 +2670,50 @@ .true. - 0 - -3.0D0 +3.22D0 3 - 1 + 2 4 - 2 + 3 1 1 - 2 - 1 + 5 4 - 4 - 7 + 2 + 4 1 2 8 + 0 + 0 + 1.0e99 1.9 1 -1 - 1.0e13 5.e15 -1 - 1.5625e13 10.e15 -1 - 1.5625e13 - 5.0e5 - 2.0e5 + 1.25e5 1.0e6 - - 0.0 - - 0.0 - 0.5 - 2 - - 0.0 - 3 - 2 0.0 -100.0 + 1.0 + + -1 + -1 + -1 1 @@ -2738,14 +2730,13 @@ 5 3 5 - 4 + 10 - 5 7 3 4 - 5 + 3 -1 -1 diff --git a/bld/namelist_files/namelist_definition.xml b/bld/namelist_files/namelist_definition.xml index ac4d3d11de..cd22e5ec0e 100644 --- a/bld/namelist_files/namelist_definition.xml +++ b/bld/namelist_files/namelist_definition.xml @@ -7210,28 +7210,6 @@ Number of hyperviscosity subcycles per dynamics timestep in sponge del2 diffusio Default: Set by build-namelist - - -Variable to specify the vertical index at which the -Rayleigh friction term is centered (the peak value). -Default: 2 - - - -Rayleigh friction parameter to determine the width of the profile. If set -to 0 then a width is chosen by the algorithm (see rayleigh_friction.F90). -Default: 0.5. - - - -Rayleigh friction parameter to determine the approximate value of the decay -time (days) at model top. If 0.0 then no Rayleigh friction is applied. -Default: 0. - - Used by SE dycore to apply sponge layer diffusion to u, v, and T for @@ -7330,6 +7308,44 @@ Second-order viscosity applied only near the model top [m^2/s]. Default: Set by build-namelist. + +Hyperviscosity coefficient se_nu [m^4/s] for u,v, T is increased to +se_nu_p*se_sponge_del4_nu_fac following a hyperbolic tangent function +centered around pressure at vertical index se_sponge_del4_lev: + + 0.5_r8*(1.0_r8+tanh(2.0_r8*log(pmid(se_sponge_del4_lev)/press))) + +where press is pressure + +If < 0, se_sponge_del4_nu_fac is automatically set based on model top location. +Default: Set by build-namelist. + + + +Divergence damping hyperviscosity coefficient se_nu_div [m^4/s] for u,v is increased to +se_nu_p*se_sponge_del4_nu_div_fac following a hyperbolic tangent function +centered around pressure at vertical index se_sponge_del4_lev: + + 0.5_r8*(1.0_r8+tanh(2.0_r8*log(pmid(se_sponge_del4_lev)/press))) + +where press is pressure + +If < 0, se_sponge_del4_nu_div_fac is automatically set based on model top location. +Default: Set by build-namelist. + + + +Level index around which increased del4 damping is centered. + +See se_sponge_del4_nu_fac and se_sponge_del4_nu_div_fac + +If < 0, se_sponge_del4_lev is automatically set based on model top location. +Default: Set by build-namelist. + + Hyperscosity for T and dp is applied to (T-Tref) and (dp-dp_ref) where diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_wcm_ne30/user_nl_cam b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_wcm_ne30/user_nl_cam index 61dcade99c..f695957589 100644 --- a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_wcm_ne30/user_nl_cam +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_wcm_ne30/user_nl_cam @@ -4,9 +4,5 @@ nhtfrq=9,9,9,9,9,9,9,9,9 inithist='ENDOFRUN' se_nsplit = 30 -se_rayk0= 10 -se_raykrange= 5 -se_raytau0=0.002 - state_debug_checks = .true. diff --git a/doc/ChangeLog b/doc/ChangeLog index 89ca2e853d..c9d150e648 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,4 +1,227 @@ =============================================================== +Tag name: cam6_3_028 +Originator(s): pel,jet +Date: Aug 19, 2021 +One-line Summary:Science Updates for CESM2.2 release and cam_development +Github PR URL:https://github.com/ESCOMP/CAM/pull/401 + +Purpose of changes (include the issue number and title text for each relevant GitHub issue): +- These mods are necessary to fix science bugs, improve algorithms and numerical stability of SE dycore. + Namelist default changes to the gravity wave parameterization and updates to the Test Tracers will + also affect the other dycores. This PR will close issue #392. + +Describe any changes made to build system: + +Describe any changes made to the namelist: +- namelist variables have been removed: + se_raytau0 + se_raykrange + se_rayk0 +- namelist variables with new defaults + effgw_beres_dp defaults + se_hypervis_scaling + se_hypervis_subcycle + se_hypervis_subcycle_sponge + se_nu + se_nu_div + se_nu_p + se_nu_top + se_molecular_diff + se_nsplit + se_rsplit + cld_macmic_num_steps +- namelist variables added for better sponge layer setup + se_sponge_del4_nu_fac + . Hyperviscosity coefficient se_nu [m^4/s] for u,v, T is increased to + se_nu_p*se_sponge_del4_nu_fac following a hyperbolic tangent function + centered around pressure at vertical index se_sponge_del4_lev: + 0.5_r8*(1.0_r8+tanh(2.0_r8*log(pmid(se_sponge_del4_lev)/press))) + where press is pressure + If less than 0, se_sponge_del4_nu_fac is automatically set based on model top location. + se_sponge_del4_nu_div_fac + . Divergence damping hyperviscosity coefficient se_nu_div [m^4/s] for u,v is increased to + se_nu_p*se_sponge_del4_nu_div_fac following a hyperbolic tangent function + centered around pressure at vertical index se_sponge_del4_lev: + 0.5_r8*(1.0_r8+tanh(2.0_r8*log(pmid(se_sponge_del4_lev)/press))) + where press is pressure + If less than 0, se_sponge_del4_nu_div_fac is automatically set based on model top location. + se_sponge_del4_lev + . Level index around which increased del4 damping is centered. + +List any changes to the defaults for the boundary datasets: + +Describe any substantial timing or memory changes: + +Code reviewed by: pel,fvitt,nusbaume,peverwhee,gold2718,cacraigucar + +List all files eliminated: + +List all files added and what they do: + +List all existing files that have been modified, and describe the changes: + +. bld/build-namelist + - Remove SE Rayleigh Friction Parameterization +. bld/namelist_files/namelist_defaults_cam.xml + - Thermodynamic consistency between pressure based physics and height based dynamics +. bld/namelist_files/namelist_definition.xml + - Remove SE Rayleigh Friction Parameterization + - New namelist variables for sponge layer configuration +. cime_config/testdefs/testmods_dirs/cam/outfrq9s_wcm_ne30/user_nl_cam + - Remove SE Rayleigh Friction Parameterization +. src/dynamics/se/dp_coupling.F90 + - Thermodynamic consistency between pressure based physics and height based dynamics +. src/dynamics/se/dycore/control_mod.F90 + - Thermodynamic consistency between pressure based physics and height based dynamics +. src/dynamics/se/dycore/dimensions_mod.F90 + - new del4 T damping variable and removed unused variables +. src/dynamics/se/dycore/global_norms_mod.F90 + - improve level dependent del4 (sponge layer) damping for hitop +. src/dynamics/se/dycore/prim_advance_mod.F90 + - Add reference temperature profile correction + - Remove SE Rayleigh Friction Parameterization + - Updated defaults for full physics, new sponge layer setup (new se_sponge... namelists parameters) +. src/dynamics/se/dycore/viscosity_mod.F90 + - add p correction to approximate Laplace on pressure surfaces +. src/dynamics/se/dyn_comp.F90 + - Remove SE Rayleigh Friction Parameterization +. src/physics/cam/check_energy.F90 + - initialize angular momentum diagnostic variables. +. src/physics/cam/geopotential.F90 + - fix dycore dependent calculation of hydrostatic elements +. src/physics/cam/physpkg.F90 + - new angular moment diagnostics +. src/physics/cam/tracers.F90 + - Changes to TT_CCOSB TT_COSB and TT_lCCOSB for idealized testing of + dynamical cores (change to Lauritzen and Thuburn,2012, QJRMS, formulation). +. src/physics/simple/physpkg.F90 + - new angular moment diagnostics +. src/utils/physconst.F90 + - correct mid-level pressure computation for SE in physconst.F90 + +If there were any failures reported from running test_driver.sh on any test +platform, and checkin with these failures has been OK'd by the gatekeeper, +then copy the lines from the td.*.status files for the failed tests to the +appropriate machine below. All failed tests must be justified. + +cheyenne/intel/aux_cam: +All tests running. Differences Expected due to namelist updates + + DIFF ERC_D_Ln9.f19_f19_mg17.QPX2000.cheyenne_intel.cam-outfrq3s + DIFF ERP_Lh12.f19_f19_mg17.FW4madSD.cheyenne_intel.cam-outfrq3h + DIFF ERP_Ln9.f19_f19_mg17.FWsc1850.cheyenne_intel.cam-outfrq9s + DIFF ERS_Ln9.f19_f19_mg17.FXSD.cheyenne_intel.cam-outfrq9s + DIFF SMS_D_Ln9.f19_f19_mg17.FWma2000climo.cheyenne_intel.cam-outfrq9s + DIFF SMS_D_Ln9.f19_f19_mg17.FXHIST.cheyenne_intel.cam-outfrq9s_amie + BASE: effgw_beres_dp = 0.1D0 + COMP: effgw_beres_dp = 0.7D0 + + DIFF ERC_D_Ln9.mpasa120z32_mpasa120.FKESSLER.cheyenne_intel.cam-outfrq3s_usecase + test tracer differences + + DIFF ERC_D_Ln9.ne16_ne16_mg17.FADIAB.cheyenne_intel.cam-terminator + DIFF ERC_D_Ln9.ne16_ne16_mg17.QPC5HIST.cheyenne_intel.cam-outfrq3s_usecase + DIFF ERP_Ln9.ne30_ne30_mg17.FCnudged.cheyenne_intel.cam-outfrq9s + DIFF SMS_Ld1.ne30pg3_ne30pg3_mg17.FC2010climo.cheyenne_intel.cam-outfrq1d + DIFF ERC_D_Ln9_P144x1.ne16pg3_ne16pg3_mg17.QPC6HIST.cheyenne_intel.cam-outfrq3s_ttrac_usecase + BASE: se_nu_top = 5.0e5 + COMP: se_nu_top = 1.25e5 + + DIFF ERP_Ln9.ne30pg3_ne30pg3_mg17.FW2000climo.cheyenne_intel.cam-outfrq9s_wcm_ne30 + BASE: se_hypervis_subcycle = 1 + COMP: se_hypervis_subcycle = 2 + BASE: se_molecular_diff = 100.0 + COMP: se_molecular_diff = 0.0 + BASE: se_nu_top = 0.0 + COMP: se_nu_top = 1.25e5 + BASE: effgw_beres_dp = 0.325D0 + COMP: effgw_beres_dp = 0.5D0 + + DIFF ERS_Ln9.ne0TESTONLYne5x4_ne0TESTONLYne5x4_mg37.FADIAB.cheyenne_intel.cam-outfrq3s_refined + BASE: se_hypervis_scaling = 3.0D0 + COMP: se_hypervis_scaling = 3.22D0 + BASE: se_hypervis_subcycle = 2 + COMP: se_hypervis_subcycle = 3 + BASE: se_hypervis_subcycle_sponge = 4 + COMP: se_hypervis_subcycle_sponge = 2 + BASE: se_nu_top = 2.0e5 + COMP: se_nu_top = 1.25e5 + + DIFF SMS_D_Ln9.ne0CONUSne30x8_ne0CONUSne30x8_mt12.FCHIST.cheyenne_intel.cam-outfrq9s_refined_camchem + BASE: se_hypervis_scaling = 3.0D0 + COMP: se_hypervis_scaling = 3.22D0 + BASE: se_hypervis_subcycle = 2 + COMP: se_hypervis_subcycle = 3 + BASE: se_hypervis_subcycle_sponge = 4 + COMP: se_hypervis_subcycle_sponge = 2 + BASE: se_nu_top = 2.0e5 + COMP: se_nu_top = 1.25e5 + BASE: cld_macmic_num_steps = 3 + COMP: cld_macmic_num_steps = 1 + + DIFF SMS_D_Ln9.ne16_ne16_mg17.QPX2000.cheyenne_intel.cam-outfrq9s + DIFF SMS_D_Ln9.ne16_ne16_mg17.FX2000.cheyenne_intel.cam-outfrq9s + BASE: se_hypervis_subcycle = 1 + COMP: se_hypervis_subcycle = 2 + BASE: se_molecular_diff = 0.0 + COMP: se_molecular_diff = 1.0 + BASE: effgw_beres_dp = 0.1D0 + COMP: effgw_beres_dp = 0.7D0 + BASE: se_hypervis_subcycle_sponge = 2 + COMP: se_hypervis_subcycle_sponge = 5 + +izumi/nag/aux_cam: all BFB except: + FAIL DAE.f45_f45_mg37.FHS94.izumi_nag.cam-dae COMPARE_base_da + - pre-existing failure + DIFF ERC_D_Ln9.f10_f10_mg37.QPWmaC6.izumi_nag.cam-outfrq3s + DIFF ERC_D_Ln9.ne16_ne16_mg17.QPC4.izumi_nag.cam-outfrq3s_usecase + DIFF ERC_D_Ln9.ne16pg3_ne16pg3_mg17.QPC4.izumi_nag.cam-outfrq3s_usecase + DIFF ERC_D_Ln9.ne5_ne5_mg37.QPC5.izumi_nag.cam-outfrq3s_ttrac + DIFF ERI_D_Ln18.ne5_ne5_mg37.FADIAB.izumi_nag.cam-outfrq3s_bwic + DIFF ERI_D_Ln18.ne5pg3_ne5pg3_mg37.FADIAB.izumi_nag.cam-outfrq3s_bwic + DIFF ERP_Ln9.ne5pg3_ne5pg3_mg37.QPC6.izumi_nag.cam-outfrq9s_clubbmf + DIFF ERS_Ln27.ne5pg3_ne5pg3_mg37.FKESSLER.izumi_nag.cam-outfrq9s + DIFF ERS_Ln9.ne5_ne5_mg37.FADIAB.izumi_nag.cam-outfrq9s + DIFF PEM_D_Ln9.ne5_ne5_mg37.FADIAB.izumi_nag.cam-outfrq3s + DIFF PLB_D_Ln9.ne5_ne5_mg37.QPC5.izumi_nag.cam-ttrac_loadbal0 + DIFF PLB_D_Ln9.ne5_ne5_mg37.QPC5.izumi_nag.cam-ttrac_loadbal1 + DIFF PLB_D_Ln9.ne5_ne5_mg37.QPC5.izumi_nag.cam-ttrac_loadbal3 + DIFF SMS_D_Ln9.ne5_ne5_mg37.QPC4X.izumi_nag.cam-outfrq9s + DIFF SMS_D_Ln9_P1x1.ne5_ne5_mg37.FADIAB.izumi_nag.cam-outfrq3s + DIFF SMS_P48x1_D_Ln9.f19_f19_mg17.FW4madSD.izumi_nag.cam-outfrq9s + - Diff due to namelist updates - see comments for cheyenne testing + +izumi/pgi/aux_cam: + DIFF ERC_D_Ln9.f10_f10_mg37.QPC4.izumi_pgi.cam-outfrq3s_diags + DIFF ERC_D_Ln9.ne5_ne5_mg37.QPC5.izumi_pgi.cam-outfrq3s_ba + DIFF ERC_D_Ln9.ne5pg2_ne5pg2_mg37.FADIAB.izumi_pgi.cam-outfrq3s + DIFF ERC_D_Ln9.ne5pg4_ne5pg4_mg37.FADIAB.izumi_pgi.cam-outfrq3s + DIFF ERP_Ln9.ne5_ne5_mg37.FHS94.izumi_pgi.cam-outfrq9s + DIFF ERP_Ln9.ne5_ne5_mg37.QPC5.izumi_pgi.cam-outfrq9s + DIFF PEM_D_Ln9.ne5pg3_ne5pg3_mg37.FADIAB.izumi_pgi.cam-outfrq3s + DIFF PLB_D_Ln9.ne5pg3_ne5pg3_mg37.QPC5.izumi_pgi.cam-ttrac_loadbal0 + DIFF PLB_D_Ln9.ne5pg3_ne5pg3_mg37.QPC5.izumi_pgi.cam-ttrac_loadbal1 + DIFF PLB_D_Ln9.ne5pg3_ne5pg3_mg37.QPC5.izumi_pgi.cam-ttrac_loadbal3 + DIFF SMS_D_Ln9.ne5pg3_ne5pg3_mg37.QPC5.izumi_pgi.cam-outfrq3s_ttrac + - Diff due to namelist updates - see comments for cheyenne testing + +Summarize any changes to answers, i.e., +- what code configurations: Configurations using SE and FV dycores +- what platforms/compilers: All +- nature of change: new climate + +If this tag changes climate describe the run(s) done to evaluate the new +climate in enough detail that it(they) could be reproduced, i.e., +- source: git clone https://github.com/jtruesdal/CAM-1.git cam6_3_019_plus_CESM2.2 + cd cam6_3_019_plus_CESM2.2 + git checkout cam6_3_019_plus_CESM2.2 +- platform/compilers: cheyenne/intel +- archived case directory: /glade/p/cesmdata/cseg/runs/cesm2_0/f.e21.FWscHIST.f09_L32_cam6_3_019_plus_CESM2.2.001.hf + +URL for AMWG diagnostics output used to validate new climate: + https://webext.cgd.ucar.edu/FWscHIST/f.e21.FWscHIST.ne30_L32_cam6_3_019_plus_CESM2.2.001.hf/atm/f.e21.FWscHIST.ne30_L32_cam6_3_019_plus_CESM2.2.001.hf_yrs1980-1989-f.e21.FWscHIST.ne30_L32_cam6_3_019.001.hf_yrs1980-1989/ +=============================================================== +=============================================================== Tag name: cam6_3_027 Originator(s): fvitt Date: 20 Jul 2021 @@ -333,7 +556,7 @@ MSS location of control simulations used to validate new climate: URL for AMWG diagnostics output used to validate new climate: - +=============================================================== =============================================================== Tag name: cam6_3_025 diff --git a/src/dynamics/se/dp_coupling.F90 b/src/dynamics/se/dp_coupling.F90 index 767c5aaaae..1b97ca1ce7 100644 --- a/src/dynamics/se/dp_coupling.F90 +++ b/src/dynamics/se/dp_coupling.F90 @@ -378,7 +378,7 @@ subroutine p_d_coupling(phys_state, phys_tend, dyn_in, tl_f, tl_qdp) end do end do call thermodynamic_consistency( & - phys_state(lchnk), phys_tend(lchnk), ncols, pver) + phys_state(lchnk), phys_tend(lchnk), ncols, pver, lchnk) end do call t_startf('pd_copy') @@ -538,7 +538,7 @@ subroutine derived_phys_dry(phys_state, phys_tend, pbuf2d) ! Finally compute energy and water column integrals of the physics input state. use constituents, only: qmin - use physconst, only: cpair, gravit, zvir, cappa, rairv, physconst_update + use physconst, only: cpairv, gravit, zvir, cappav, rairv, physconst_update use shr_const_mod, only: shr_const_rwv use phys_control, only: waccmx_is use geopotential, only: geopotential_t @@ -644,8 +644,6 @@ subroutine derived_phys_dry(phys_state, phys_tend, pbuf2d) do k = 1, nlev do i = 1, ncol phys_state(lchnk)%rpdel(i,k) = 1._r8/phys_state(lchnk)%pdel(i,k) - phys_state(lchnk)%exner (i,k) = (phys_state(lchnk)%pint(i,pver+1) & - / phys_state(lchnk)%pmid(i,k))**cappa end do end do @@ -701,12 +699,20 @@ subroutine derived_phys_dry(phys_state, phys_tend, pbuf2d) ! Fill local zvirv variable; calculated for WACCM-X. !----------------------------------------------------------------------------- if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then - call physconst_update(phys_state(lchnk)%q, phys_state(lchnk)%t, lchnk, ncol) + call physconst_update(phys_state(lchnk)%q, phys_state(lchnk)%t, lchnk, ncol,& + to_moist_factor=phys_state(lchnk)%pdeldry(:ncol,:)/phys_state(lchnk)%pdel(:ncol,:) ) zvirv(:,:) = shr_const_rwv / rairv(:,:,lchnk) -1._r8 else zvirv(:,:) = zvir endif + do k = 1, nlev + do i = 1, ncol + phys_state(lchnk)%exner(i,k) = (phys_state(lchnk)%pint(i,pver+1) & + / phys_state(lchnk)%pmid(i,k))**cappav(i,k,lchnk) + end do + end do + ! Compute initial geopotential heights - based on full pressure call geopotential_t (phys_state(lchnk)%lnpint, phys_state(lchnk)%lnpmid , phys_state(lchnk)%pint , & phys_state(lchnk)%pmid , phys_state(lchnk)%pdel , phys_state(lchnk)%rpdel , & @@ -716,8 +722,8 @@ subroutine derived_phys_dry(phys_state, phys_tend, pbuf2d) ! Compute initial dry static energy, include surface geopotential do k = 1, pver do i = 1, ncol - phys_state(lchnk)%s(i,k) = cpair*phys_state(lchnk)%t(i,k) & - + gravit*phys_state(lchnk)%zm(i,k) + phys_state(lchnk)%phis(i) + phys_state(lchnk)%s(i,k) = cpairv(i,k,lchnk)*phys_state(lchnk)%t(i,k) & + + gravit*phys_state(lchnk)%zm(i,k) + phys_state(lchnk)%phis(i) end do end do @@ -736,19 +742,20 @@ end subroutine derived_phys_dry !========================================================================================= -subroutine thermodynamic_consistency(phys_state, phys_tend, ncols, pver) +subroutine thermodynamic_consistency(phys_state, phys_tend, ncols, pver, lchnk) ! ! Adjust the physics temperature tendency for thermal energy consistency with the ! dynamics. ! Note: mixing ratios are assumed to be dry. ! use dimensions_mod, only: lcp_moist - use physconst, only: get_cp, cpair + use physconst, only: get_cp use control_mod, only: phys_dyn_cp + use physconst, only: cpairv type(physics_state), intent(in) :: phys_state type(physics_tend ), intent(inout) :: phys_tend - integer, intent(in) :: ncols, pver + integer, intent(in) :: ncols, pver, lchnk real(r8):: inv_cp(ncols,pver) !---------------------------------------------------------------------------- @@ -762,7 +769,7 @@ subroutine thermodynamic_consistency(phys_state, phys_tend, ncols, pver) ! consistency (not taking into account dme adjust) ! call get_cp(1,ncols,1,pver,1,1,pcnst,phys_state%q(1:ncols,1:pver,:),.true.,inv_cp) - phys_tend%dtdt(1:ncols,1:pver) = phys_tend%dtdt(1:ncols,1:pver)*cpair*inv_cp + phys_tend%dtdt(1:ncols,1:pver) = phys_tend%dtdt(1:ncols,1:pver)*cpairv(1:ncols,1:pver,lchnk)*inv_cp end if end subroutine thermodynamic_consistency diff --git a/src/dynamics/se/dycore/control_mod.F90 b/src/dynamics/se/dycore/control_mod.F90 index d5fc4abe81..0ecc2079d5 100644 --- a/src/dynamics/se/dycore/control_mod.F90 +++ b/src/dynamics/se/dycore/control_mod.F90 @@ -23,7 +23,7 @@ module control_mod ! every rsplit tracer timesteps logical, public :: variable_nsplit=.false. - integer, public :: phys_dyn_cp = 0 !=0; no thermal energy scaling of T increment + integer, public :: phys_dyn_cp = 1 !=0; no thermal energy scaling of T increment !=1; scale increment for cp consistency between dynamics and physics logical, public :: refined_mesh @@ -63,10 +63,25 @@ module control_mod ! (only used for variable viscosity, recommend 1.9 in namelist) real (kind=r8), public :: nu = 7.0D5 ! viscosity (momentum equ) real (kind=r8), public :: nu_div = -1 ! viscsoity (momentum equ, div component) - real (kind=r8), public :: nu_s = -1 ! default = nu T equ. viscosity + real (kind=r8), public :: nu_t = -1 ! default = nu T equ. viscosity real (kind=r8), public :: nu_q = -1 ! default = nu tracer viscosity real (kind=r8), public :: nu_p = 0.0D5 ! default = 0 ps equ. viscosity real (kind=r8), public :: nu_top = 0.0D5 ! top-of-the-model viscosity + + ! + ! Del4 sponge layer diffusion + ! + ! Divergence damping hyperviscosity coefficient nu_div [m^4/s] for u,v is increased to + ! nu_div*sponge_del4_nu_div_fac following a hyperbolic tangent function + ! centered around pressure at vertical index sponge_del4_lev + ! + ! Similar for sponge_del4_nu_fac + ! + real(r8), public :: sponge_del4_nu_fac + real(r8), public :: sponge_del4_nu_div_fac + integer , public :: sponge_del4_lev + + integer, public :: hypervis_subcycle=1 ! number of subcycles for hyper viscsosity timestep integer, public :: hypervis_subcycle_sponge=1 ! number of subcycles for hyper viscsosity timestep in sponge integer, public :: hypervis_subcycle_q=1 ! number of subcycles for hyper viscsosity timestep on TRACERS @@ -106,12 +121,6 @@ module control_mod integer, public, parameter :: nwest = 7 integer, public, parameter :: neast = 8 - ! - ! parameters for sponge layer Rayleigh damping - ! - real(r8), public :: raytau0 - real(r8), public :: raykrange - integer, public :: rayk0 ! ! molecular diffusion ! diff --git a/src/dynamics/se/dycore/dimensions_mod.F90 b/src/dynamics/se/dycore/dimensions_mod.F90 index a012c76148..5d3d52a448 100644 --- a/src/dynamics/se/dycore/dimensions_mod.F90 +++ b/src/dynamics/se/dycore/dimensions_mod.F90 @@ -73,8 +73,8 @@ module dimensions_mod integer, allocatable, public :: kord_tr(:), kord_tr_cslam(:) real(r8), public :: nu_scale_top(PLEV)! scaling of del2 viscosity in sopnge layer (initialized in dyn_comp) - real(r8), public :: nu_lev(PLEV) - real(r8), public :: otau(PLEV) + real(r8), public :: nu_lev(PLEV) ! level dependent del4 (u,v) damping + real(r8), public :: nu_t_lev(PLEV) ! level depedendet del4 T damping integer, public :: ksponge_end ! sponge is active k=1,ksponge_end real(r8), public :: nu_div_lev(PLEV) = 1.0_r8 ! scaling of viscosity in sponge layer ! (set in prim_state; if applicable) @@ -84,7 +84,6 @@ module dimensions_mod real(r8), public :: km_sponge_factor(PLEV) !scaling for molecular diffusion (when used as sponge) real(r8), public :: kmvisi_ref(PLEV+1) !reference profiles for molecular diffusion real(r8), public :: kmcndi_ref(PLEV+1) !reference profiles for molecular diffusion - real(r8), public :: rhoi_ref(PLEV+1) !reference profiles for rho integer, public :: nhc_phys diff --git a/src/dynamics/se/dycore/global_norms_mod.F90 b/src/dynamics/se/dycore/global_norms_mod.F90 index 8f55f6391e..49f104380e 100644 --- a/src/dynamics/se/dycore/global_norms_mod.F90 +++ b/src/dynamics/se/dycore/global_norms_mod.F90 @@ -208,14 +208,15 @@ subroutine print_cfl(elem,hybrid,nets,nete,dtnu,ptop,pmid,& use hybrid_mod, only: hybrid_t, PrintHybrid use element_mod, only: element_t use dimensions_mod, only: np,ne,nelem,nelemd,nc,nhe,qsize,ntrac,nlev,large_Courant_incr - use dimensions_mod, only: nu_scale_top,nu_div_lev,nu_lev + use dimensions_mod, only: nu_scale_top,nu_div_lev,nu_lev,nu_t_lev use quadrature_mod, only: gausslobatto, quadrature_t use reduction_mod, only: ParallelMin,ParallelMax use physconst, only: ra, rearth, pi - use control_mod, only: nu, nu_div, nu_q, nu_p, nu_s, nu_top, fine_ne, rk_stage_user, max_hypervis_courant + use control_mod, only: nu, nu_div, nu_q, nu_p, nu_t, nu_top, fine_ne, rk_stage_user, max_hypervis_courant use control_mod, only: tstep_type, hypervis_power, hypervis_scaling + use control_mod, only: sponge_del4_nu_div_fac, sponge_del4_nu_fac, sponge_del4_lev use cam_abortutils, only: endrun use parallel_mod, only: global_shared_buf, global_shared_sum use edge_mod, only: initedgebuffer, FreeEdgeBuffer, edgeVpack, edgeVunpack @@ -246,7 +247,7 @@ subroutine print_cfl(elem,hybrid,nets,nete,dtnu,ptop,pmid,& real (kind=r8) :: x, y, noreast, nw, se, sw real (kind=r8), dimension(np,np,nets:nete) :: zeta real (kind=r8) :: lambda_max, lambda_vis, min_gw, lambda,umax, ugw - real (kind=r8) :: press,scale1,scale2,scale3, max_laplace + real (kind=r8) :: scale1,scale2,scale3, max_laplace integer :: ie,corner, i, j, rowind, colind, k type (quadrature_t) :: gp character(LEN=256) :: rk_str @@ -255,11 +256,12 @@ subroutine print_cfl(elem,hybrid,nets,nete,dtnu,ptop,pmid,& real (kind=r8) :: dt_max_adv, dt_max_gw, dt_max_tracer_se, dt_max_tracer_fvm real (kind=r8) :: dt_max_hypervis, dt_max_hypervis_tracer, dt_max_laplacian_top - real(kind=r8) :: I_sphere + real(kind=r8) :: I_sphere, nu_max, nu_div_max real(kind=r8) :: h(np,np,nets:nete) + logical :: top_000_042km, top_042_090km, top_090_140km, top_140_600km ! model top location ranges + logical :: nu_set,div_set,lev_set - ! Eigenvalues calculated by folks at UMich (Paul U & Jared W) select case (np) case (2) @@ -549,34 +551,100 @@ subroutine print_cfl(elem,hybrid,nets,nete,dtnu,ptop,pmid,& call automatically_set_viscosity_coefficients(hybrid,ne,max_min_dx,min_min_dx,nu_p ,1.0_r8 ,'_p ') call automatically_set_viscosity_coefficients(hybrid,ne,max_min_dx,min_min_dx,nu ,0.5_r8,' ') - if (ptop>100.0_r8) then + call automatically_set_viscosity_coefficients(hybrid,ne,max_min_dx,min_min_dx,nu_div,2.5_r8 ,'_div') + + if (nu_q<0) nu_q = nu_p ! necessary for consistency + if (nu_t<0) nu_t = nu_p ! temperature damping is always equal to nu_p + + nu_div_lev(:) = nu_div + nu_lev(:) = nu + nu_t_lev(:) = nu_p + + ! + ! sponge layer strength needed for stability depends on model top location + ! + top_000_042km = .false. + top_042_090km = .false. + top_090_140km = .false. + top_140_600km = .false. + nu_set = sponge_del4_nu_fac < 0 + div_set = sponge_del4_nu_div_fac < 0 + lev_set = sponge_del4_lev < 0 + if (ptop>1.0_r8) then + ! + ! CAM6 top (~2.3 Pa) ! - ! CAM setting + top_000_042km = .true. + else if (ptop>3e-3_r8) then ! - call automatically_set_viscosity_coefficients(hybrid,ne,max_min_dx,min_min_dx,nu_div,2.5_r8 ,'_div') - nu_div_lev(:) = nu_div - nu_lev(:) = nu + ! CAM7 top (~4.35e-3 Pa) + ! + top_042_090km = .true. + else if (ptop>3E-6_r8) then + ! + ! WACCM top (~4.5e-6 Pa) + ! + top_090_140km = .true. else ! - ! WACCM setting + ! WACCM-x - geospace ! - call automatically_set_viscosity_coefficients(hybrid,ne,max_min_dx,min_min_dx,nu_div,2.5_r8 ,'_div') - if (hybrid%masterthread) write(iulog,*) ": sponge layer viscosity scaling factor" - do k=1,nlev - press = pmid(k) - - scale1 = 0.5_r8*(1.0_r8+tanh(2.0_r8*log(100.0_r8/press))) - nu_div_lev(k) = (1.0_r8-scale1)*nu_div+scale1*2.0_r8*nu_div - nu_div_lev(k) = nu_div - nu_lev(k) = (1.0_r8-scale1)*nu +scale1*nu_p - nu_lev(k) = nu - if (hybrid%masterthread) write(iulog,*) "nu_lev=",k,nu_lev(k) - if (hybrid%masterthread) write(iulog,*) "nu_div_lev=",k,nu_div_lev(k) - end do + top_140_600km = .true. + end if + ! + ! Logging text for sponge layer configuration + ! + if (hybrid%masterthread .and. nu_set .or. div_set .or. lev_set) then + write(iulog,* )"" + write(iulog,* )"Sponge layer del4 coefficient defaults based on model top location:" + end if + ! + ! if user or namelist is not specifying sponge del4 settings here are best guesses (empirically determined) + ! + if (top_000_042km) then + if (sponge_del4_lev <0) sponge_del4_lev = 3 + if (sponge_del4_nu_fac <0) sponge_del4_nu_fac = 1.0_r8 + if (sponge_del4_nu_div_fac<0) sponge_del4_nu_div_fac = 4.5_r8 end if - if (nu_q<0) nu_q = nu_p ! necessary for consistency - if (nu_s<0) nu_s = nu_p ! temperature damping is always equal to nu_p + if (top_042_090km) then + if (sponge_del4_lev <0) sponge_del4_lev = 3 + if (sponge_del4_nu_fac <0) sponge_del4_nu_fac = 5.0_r8 + if (sponge_del4_nu_div_fac<0) sponge_del4_nu_div_fac = 7.5_r8 + end if + + if (top_090_140km.or.top_140_600km) then + if (sponge_del4_lev <0) sponge_del4_lev = 10 + if (sponge_del4_nu_fac <0) sponge_del4_nu_fac = 5.0_r8 + if (sponge_del4_nu_div_fac<0) sponge_del4_nu_div_fac = 7.5_r8 + end if + ! + ! Log sponge layer configuration + ! + if (hybrid%masterthread.and.nu_set) & + write(iulog, '(a,e9.2)') ' sponge_del4_nu_fac = ',sponge_del4_nu_fac + if (hybrid%masterthread.and.div_set) & + write(iulog, '(a,e9.2)') ' sponge_del4_nu_div_fac = ',sponge_del4_nu_div_fac + if (hybrid%masterthread.and.lev_set) & + write(iulog, '(a,i0)') ' sponge_del4_lev = ',sponge_del4_lev + write(iulog,* )"" + + if (hybrid%masterthread) write(iulog,*) ": sponge layer viscosity scaling factor" + nu_max = sponge_del4_nu_fac*nu_p + nu_div_max = sponge_del4_nu_div_fac*nu_p + do k=1,nlev + ! Vertical profile from FV dycore (see Lauritzen et al. 2012 DOI:10.1177/1094342011410088) + scale1 = 0.5_r8*(1.0_r8+tanh(2.0_r8*log(pmid(sponge_del4_lev)/pmid(k)))) + nu_div_lev(k) = (1.0_r8-scale1)*nu_div+scale1*nu_div_max + if (sponge_del4_nu_fac.ne.1.0_r8) then + nu_lev(k) = (1.0_r8-scale1)*nu +scale1*nu_max + nu_t_lev(k) = (1.0_r8-scale1)*nu_p +scale1*nu_max + end if + + if (hybrid%masterthread) write(iulog,*) "nu_t_lev =",k,nu_t_lev(k) + if (hybrid%masterthread) write(iulog,*) "nu,nu_div_lev=",k,nu_lev(k),nu_div_lev(k) + end do + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! @@ -607,9 +675,6 @@ subroutine print_cfl(elem,hybrid,nets,nete,dtnu,ptop,pmid,& write(iulog,*) rk_str write(iulog,'(a)') ' * Spectral-element advection uses SSP preservation RK3' write(iulog,'(a)') ' * Viscosity operators use forward Euler' - if (ntrac>0) then - write(iulog,'(a)') ' * CSLAM uses two time-levels backward trajectory method' - end if end if S_laplacian = 2.0_r8 !using forward Euler for sponge diffusion S_hypervis = 2.0_r8 !using forward Euler for hyperviscosity @@ -636,7 +701,8 @@ subroutine print_cfl(elem,hybrid,nets,nete,dtnu,ptop,pmid,& else dt_max_tracer_fvm = -1.0_r8 end if - dt_max_hypervis = s_hypervis/(MAX(MAXVAL(nu_div_lev(:)),MAXVAL(nu_lev(:)))*normDinv_hypervis) + nu_max = MAX(MAXVAL(nu_div_lev(:)),MAXVAL(nu_lev(:)),MAXVAL(nu_t_lev(:))) + dt_max_hypervis = s_hypervis/(nu_max*normDinv_hypervis) dt_max_hypervis_tracer = s_hypervis/(nu_q*normDinv_hypervis) max_laplace = MAX(MAXVAL(nu_scale_top(:))*nu_top,MAXVAL(kmvis_ref(:)/rho_ref(:))) @@ -676,8 +742,14 @@ subroutine print_cfl(elem,hybrid,nets,nete,dtnu,ptop,pmid,& write(iulog,'(a,f10.2,a,f10.2,a)') '* dt (del2 sponge ; u,v,T,dM) < ',& dt_max_laplacian_top,'s',dt_dyn_del2_actual,'s' - if (dt_dyn_del2_actual>dt_max_laplacian_top) & - write(iulog,*) 'WARNING: theoretically unstable in sponge; increase se_hypervis_subcycle_sponge' + if (dt_dyn_del2_actual>dt_max_laplacian_top) then + if (k==1) then + write(iulog,*) 'WARNING: theoretically unstable in sponge; increase se_hypervis_subcycle_sponge',& + ' (this WARNING can sometimes be ignored in level 1)' + else + write(iulog,*) 'WARNING: theoretically unstable in sponge; increase se_hypervis_subcycle_sponge' + endif + end if end do write(iulog,*) ' ' if (hypervis_power /= 0) then @@ -1107,7 +1179,9 @@ subroutine automatically_set_viscosity_coefficients(hybrid,ne,max_min_dx,min_min if (nu < 0) then if (ne <= 0) then - if (hypervis_scaling/=0) then + if (hypervis_power/=0) then + call endrun('ERROR: Automatic scaling of scalar viscosity not implemented') + else if (hypervis_scaling/=0) then nu_min = factor*nu_fac*(max_min_dx*1000.0_r8)**uniform_res_hypervis_scaling nu_max = factor*nu_fac*(min_min_dx*1000.0_r8)**uniform_res_hypervis_scaling nu = factor*nu_min @@ -1119,8 +1193,6 @@ subroutine automatically_set_viscosity_coefficients(hybrid,ne,max_min_dx,min_min nu = nu_min*(2.0_r8*rearth/(3.0_r8*max_min_dx*1000.0_r8))**hypervis_scaling/(rearth**4) if (hybrid%masterthread) & write(iulog,'(a,a,a,e9.3)') "Nu_tensor",TRIM(str)," = ",nu - else if (hypervis_power/=0) then - call endrun('ERROR: Automatic scaling of scalar viscosity not implemented') end if else nu = factor*nu_fac*((30.0_r8/ne)*110000.0_r8)**uniform_res_hypervis_scaling diff --git a/src/dynamics/se/dycore/prim_advance_mod.F90 b/src/dynamics/se/dycore/prim_advance_mod.F90 index 2f06fe7119..459658200f 100644 --- a/src/dynamics/se/dycore/prim_advance_mod.F90 +++ b/src/dynamics/se/dycore/prim_advance_mod.F90 @@ -56,7 +56,6 @@ subroutine prim_advance_exp(elem, fvm, deriv, hvcoord, hybrid,dt, tl, nets, net use time_mod, only: TimeLevel_t, timelevel_qdp, tevolve use dimensions_mod, only: lcp_moist use fvm_control_volume_mod, only: fvm_struct - use control_mod, only: raytau0 use physconst, only: get_cp, thermodynamic_active_species_num use physconst, only: get_kappa_dry, dry_air_species_num use physconst, only: thermodynamic_active_species_idx_dycore @@ -147,7 +146,6 @@ subroutine prim_advance_exp(elem, fvm, deriv, hvcoord, hybrid,dt, tl, nets, net dt_vis = dt - if (raytau0>0) call rayleigh_friction(elem,n0,nets,nete,dt) if (tstep_type==1) then ! RK2-SSP 3 stage. matches tracer scheme. optimal SSP CFL, but ! not optimal for regular CFL @@ -441,7 +439,7 @@ subroutine advance_hypervis_dp(edge3,elem,fvm,hybrid,deriv,nt,qn0,nets,nete,dt2, ! ! take one timestep of: ! u(:,:,:,np) = u(:,:,:,np) + dt2*nu*laplacian**order ( u ) - ! T(:,:,:,np) = T(:,:,:,np) + dt2*nu_s*laplacian**order ( T ) + ! T(:,:,:,np) = T(:,:,:,np) + dt2*nu_t*laplacian**order ( T ) ! ! ! For correct scaling, dt2 should be the same 'dt2' used in the leapfrog advace @@ -451,8 +449,8 @@ subroutine advance_hypervis_dp(edge3,elem,fvm,hybrid,deriv,nt,qn0,nets,nete,dt2, use dimensions_mod, only: np, nlev, nc, ntrac, npsq, qsize use dimensions_mod, only: hypervis_dynamic_ref_state,ksponge_end use dimensions_mod, only: nu_scale_top,nu_lev,kmvis_ref,kmcnd_ref,rho_ref,km_sponge_factor - use dimensions_mod, only: kmvisi_ref,kmcndi_ref,rhoi_ref - use control_mod, only: nu, nu_s, hypervis_subcycle,hypervis_subcycle_sponge, nu_p, nu_top + use dimensions_mod, only: kmvisi_ref,kmcndi_ref,nu_t_lev + use control_mod, only: nu, nu_t, hypervis_subcycle,hypervis_subcycle_sponge, nu_p, nu_top use control_mod, only: molecular_diff use hybrid_mod, only: hybrid_t!, get_loop_ranges use element_mod, only: element_t @@ -483,7 +481,7 @@ subroutine advance_hypervis_dp(edge3,elem,fvm,hybrid,deriv,nt,qn0,nets,nete,dt2, integer :: kbeg, kend, kblk real (kind=r8), dimension(np,np,2,nlev,nets:nete) :: vtens real (kind=r8), dimension(np,np,nlev,nets:nete) :: ttens, dptens - real (kind=r8), dimension(np,np,nlev,nets:nete) :: dp3d_ref, T_ref + real (kind=r8), dimension(np,np,nlev,nets:nete) :: dp3d_ref, T_ref, pmid_ref real (kind=r8), dimension(np,np,nets:nete) :: ps_ref real (kind=r8), dimension(0:np+1,0:np+1,nlev) :: corners real (kind=r8), dimension(2,2,2) :: cflux @@ -496,8 +494,6 @@ subroutine advance_hypervis_dp(edge3,elem,fvm,hybrid,deriv,nt,qn0,nets,nete,dt2, real (kind=r8), dimension(np,np) :: tmp, tmp2 real (kind=r8), dimension(np,np,ksponge_end,nets:nete):: kmvis,kmcnd,rho_dry real (kind=r8), dimension(np,np,ksponge_end+1):: kmvisi,kmcndi - real (kind=r8), dimension(np,np,ksponge_end+1):: pint,rhoi_dry - real (kind=r8), dimension(np,np,ksponge_end ):: pmid real (kind=r8), dimension(np,np,nlev) :: tmp_kmvis,tmp_kmcnd real (kind=r8), dimension(np,np,2) :: lap_v real (kind=r8) :: v1,v2,v1new,v2new,dt,heating,T0,T1 @@ -507,7 +503,7 @@ subroutine advance_hypervis_dp(edge3,elem,fvm,hybrid,deriv,nt,qn0,nets,nete,dt2, real (kind=r8), dimension(ksponge_end) :: dtemp,du,dv real (kind=r8) :: nu_temp, nu_dp, nu_velo - if (nu_s == 0 .and. nu == 0 .and. nu_p==0 ) return; + if (nu_t == 0 .and. nu == 0 .and. nu_p==0 ) return; ptop = hvcoord%hyai(1)*hvcoord%ps0 @@ -539,11 +535,15 @@ subroutine advance_hypervis_dp(edge3,elem,fvm,hybrid,deriv,nt,qn0,nets,nete,dt2, T0 = Tref-T1 do ie=nets,nete do k=1,nlev + pmid_ref(:,:,k,ie) =hvcoord%hyam(k)*hvcoord%ps0 + hvcoord%hybm(k)*ps_ref(:,:,ie) dp3d_ref(:,:,k,ie) = ((hvcoord%hyai(k+1)-hvcoord%hyai(k))*hvcoord%ps0 + & (hvcoord%hybi(k+1)-hvcoord%hybi(k))*ps_ref(:,:,ie)) - tmp = hvcoord%hyam(k)*hvcoord%ps0+hvcoord%hybm(k)*ps_ref(:,:,ie) - tmp2 = (tmp/hvcoord%ps0)**cappa - T_ref(:,:,k,ie) = (T0+T1*tmp2) + if (hvcoord%hybm(k)>0) then + tmp2 = (pmid_ref(:,:,k,ie)/hvcoord%ps0)**cappa + T_ref(:,:,k,ie) = (T0+T1*tmp2) + else + T_ref(:,:,k,ie) = 0.0_r8 + end if end do end do @@ -561,8 +561,8 @@ subroutine advance_hypervis_dp(edge3,elem,fvm,hybrid,deriv,nt,qn0,nets,nete,dt2, call calc_tot_energy_dynamics(elem,fvm,nets,nete,nt,qn0,'dBH') rhypervis_subcycle=1.0_r8/real(hypervis_subcycle,kind=r8) - call biharmonic_wk_dp3d(elem,dptens,dpflux,ttens,vtens,deriv,edge3,hybrid,nt,nets,nete,kbeg,kend,& - dp3d_ref,T_ref) + call biharmonic_wk_dp3d(elem,dptens,dpflux,ttens,vtens,deriv,edge3,hybrid,nt,nets,nete,kbeg,kend,hvcoord,& + dp3d_ref=dp3d_ref,pmid_ref=pmid_ref) do ie=nets,nete ! compute mean flux @@ -590,7 +590,7 @@ subroutine advance_hypervis_dp(edge3,elem,fvm,hybrid,deriv,nt,qn0,nets,nete,dt2, !DIR_VECTOR_ALIGNED do j=1,np do i=1,np - ttens(i,j,k,ie) = -nu_s*ttens(i,j,k,ie) + ttens(i,j,k,ie) = -nu_t_lev(k)*ttens(i,j,k,ie) dptens(i,j,k,ie) = -nu_p*dptens(i,j,k,ie) vtens(i,j,1,k,ie) = -nu_lev(k)*vtens(i,j,1,k,ie) vtens(i,j,2,k,ie) = -nu_lev(k)*vtens(i,j,2,k,ie) @@ -749,57 +749,11 @@ subroutine advance_hypervis_dp(edge3,elem,fvm,hybrid,deriv,nt,qn0,nets,nete,dt2, ! !*************************************************************** ! - ! - ! vertical diffusion - ! - call t_startf('vertical_molec_diff') - if (molecular_diff>1) then - do ie=nets,nete - call get_rho_dry(1,np,1,np,ksponge_end,nlev,qsize,elem(ie)%state%Qdp(:,:,:,1:qsize,qn0), & - elem(ie)%state%T(:,:,:,nt),ptop,elem(ie)%state%dp3d(:,:,:,nt),& - .true.,rhoi_dry=rhoi_dry(:,:,:), & - active_species_idx_dycore=thermodynamic_active_species_idx_dycore,& - pint_out=pint,pmid_out=pmid) - ! - ! constant coefficients - ! - do k=1,ksponge_end+1 - kmvisi(:,:,k) = kmvisi_ref(k)*rhoi_dry(:,:,k) - kmcndi(:,:,k) = kmcndi_ref(k)*rhoi_dry(:,:,k) - end do - ! - ! do vertical diffusion - ! - do j=1,np - do i=1,np - call solve_diffusion(dt2,np,nlev,i,j,ksponge_end,pmid,pint,kmcndi(:,:,:)/cpair,elem(ie)%state%T(:,:,:,nt),& - 0,dtemp) - call solve_diffusion(dt2,np,nlev,i,j,ksponge_end,pmid,pint,kmvisi(:,:,:),elem(ie)%state%v(:,:,1,:,nt),1,du) - call solve_diffusion(dt2,np,nlev,i,j,ksponge_end,pmid,pint,kmvisi(:,:,:),elem(ie)%state%v(:,:,2,:,nt),1,dv) - do k=1,ksponge_end - v1 = elem(ie)%state%v(i,j,1,k,nt) - v2 = elem(ie)%state%v(i,j,2,k,nt) - v1new = v1 + du(k) - v2new = v2 + dv(k) - ! - ! frictional heating - ! - heating = 0.5_r8*((v1new*v1new+v2new*v2new) - (v1*v1+v2*v2)) - elem(ie)%state%T(i,j,k,nt)=elem(ie)%state%T(i,j,k,nt) & - -heating*inv_cp_full(i,j,k,ie)+dtemp(k) - elem(ie)%state%v(i,j,1,k,nt)=v1new - elem(ie)%state%v(i,j,2,k,nt)=v2new - end do - end do - end do - end do - end if - call t_stopf('vertical_molec_diff') call t_startf('sponge_diff') ! ! compute coefficients for horizontal diffusion ! - if (molecular_diff>0) then + if (molecular_diff==1) then do ie=nets,nete call get_rho_dry(1,np,1,np,ksponge_end,nlev,qsize,elem(ie)%state%Qdp(:,:,:,1:qsize,qn0), & elem(ie)%state%T(:,:,:,nt),ptop,elem(ie)%state%dp3d(:,:,:,nt),& @@ -807,27 +761,15 @@ subroutine advance_hypervis_dp(edge3,elem,fvm,hybrid,deriv,nt,qn0,nets,nete,dt2, active_species_idx_dycore=thermodynamic_active_species_idx_dycore) end do - if (molecular_diff==1) then - do ie=nets,nete - ! - ! compute molecular diffusion and thermal conductivity coefficients at mid-levels - ! - call get_molecular_diff_coef(1,np,1,np,ksponge_end,nlev,& - elem(ie)%state%T(:,:,:,nt),0,km_sponge_factor(1:ksponge_end),kmvis(:,:,:,ie),kmcnd(:,:,:,ie),qsize,& - elem(ie)%state%Qdp(:,:,:,1:qsize,qn0),fact=1.0_r8/elem(ie)%state%dp3d(:,:,1:ksponge_end,nt),& - active_species_idx_dycore=thermodynamic_active_species_idx_dycore) - end do - else + do ie=nets,nete ! - ! constant coefficients + ! compute molecular diffusion and thermal conductivity coefficients at mid-levels ! - do ie=nets,nete - do k=1,ksponge_end - kmvis (:,:,k,ie) = kmvis_ref(k) - kmcnd (:,:,k,ie) = kmcnd_ref(k) - end do - end do - end if + call get_molecular_diff_coef(1,np,1,np,ksponge_end,nlev,& + elem(ie)%state%T(:,:,:,nt),0,km_sponge_factor(1:ksponge_end),kmvis(:,:,:,ie),kmcnd(:,:,:,ie),qsize,& + elem(ie)%state%Qdp(:,:,:,1:qsize,qn0),fact=1.0_r8/elem(ie)%state%dp3d(:,:,1:ksponge_end,nt),& + active_species_idx_dycore=thermodynamic_active_species_idx_dycore) + end do ! ! diagnostics ! @@ -1097,7 +1039,7 @@ subroutine compute_and_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,& use physconst, only: epsilo, get_gz_given_dp_Tv_Rdry use physconst, only: thermodynamic_active_species_num, get_virtual_temp, get_cp_dry use physconst, only: thermodynamic_active_species_idx_dycore,get_R_dry - use physconst, only: dry_air_species_num,get_exner + use physconst, only: dry_air_species_num,get_exner,tref,cpair,gravit,lapse_rate use time_mod, only : tevolve implicit none @@ -1144,6 +1086,9 @@ subroutine compute_and_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,& real (kind=r8) :: stashdp3d (np,np,nlev),tempdp3d(np,np), tempflux(nc,nc,4) real (kind=r8) :: ckk, term, T_v(np,np,nlev) real (kind=r8), dimension(np,np,2) :: grad_exner + real (kind=r8), dimension(np,np,2) :: grad_exner_term + real (kind=r8), dimension(np,np,2) :: grad_logexner + real (kind=r8) :: T0,T1 real (kind=r8), dimension(np,np) :: theta_v type (EdgeDescriptor_t):: desc @@ -1295,8 +1240,8 @@ subroutine compute_and_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,& theta_v(:,:)=T_v(:,:,k)/exner(:,:) call gradient_sphere(exner(:,:),deriv,elem(ie)%Dinv,grad_exner) - grad_exner(:,:,1) = cp_dry(:,:,k)*theta_v(:,:)*grad_exner(:,:,1) - grad_exner(:,:,2) = cp_dry(:,:,k)*theta_v(:,:)*grad_exner(:,:,2) + grad_exner_term(:,:,1) = cp_dry(:,:,k)*theta_v(:,:)*grad_exner(:,:,1) + grad_exner_term(:,:,2) = cp_dry(:,:,k)*theta_v(:,:)*grad_exner(:,:,2) else exner(:,:)=(p_full(:,:,k)/hvcoord%ps0)**kappa(:,:,k,ie) theta_v(:,:)=T_v(:,:,k)/exner(:,:) @@ -1307,14 +1252,34 @@ subroutine compute_and_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,& grad_kappa_term(:,:,1)=-suml*grad_kappa_term(:,:,1) grad_kappa_term(:,:,2)=-suml*grad_kappa_term(:,:,2) - grad_exner(:,:,1) = cp_dry(:,:,k)*theta_v(:,:)*(grad_exner(:,:,1)+grad_kappa_term(:,:,1)) - grad_exner(:,:,2) = cp_dry(:,:,k)*theta_v(:,:)*(grad_exner(:,:,2)+grad_kappa_term(:,:,2)) + grad_exner_term(:,:,1) = cp_dry(:,:,k)*theta_v(:,:)*(grad_exner(:,:,1)+grad_kappa_term(:,:,1)) + grad_exner_term(:,:,2) = cp_dry(:,:,k)*theta_v(:,:)*(grad_exner(:,:,2)+grad_kappa_term(:,:,2)) + end if + + ! balanced ref profile correction: + ! reference temperature profile (Simmons and Jiabin, 1991, QJRMS, Section 2a) + ! + ! Tref = T0+T1*Exner + ! T1 = .0065*Tref*Cp/g ! = ~191 + ! T0 = Tref-T1 ! = ~97 + ! + T1 = lapse_rate*Tref*cpair/gravit + T0 = Tref-T1 + + if (hvcoord%hybm(k)>0) then + call gradient_sphere(log(exner(:,:)),deriv,elem(ie)%Dinv,grad_logexner) + grad_exner_term(:,:,1)=grad_exner_term(:,:,1) + & + cpair*T0*(grad_logexner(:,:,1)-grad_exner(:,:,1)/exner(:,:)) + grad_exner_term(:,:,2)=grad_exner_term(:,:,2) + & + cpair*T0*(grad_logexner(:,:,2)-grad_exner(:,:,2)/exner(:,:)) end if + + do j=1,np do i=1,np - glnps1 = grad_exner(i,j,1) - glnps2 = grad_exner(i,j,2) + glnps1 = grad_exner_term(i,j,1) + glnps2 = grad_exner_term(i,j,2) v1 = elem(ie)%state%v(i,j,1,k,n0) v2 = elem(ie)%state%v(i,j,2,k,n0) @@ -2179,78 +2144,4 @@ subroutine fill_element(Eval) return end subroutine fill_element - subroutine rayleigh_friction(elem,nt,nets,nete,dt) - use dimensions_mod, only: nlev, otau - use hybrid_mod, only: hybrid_t!, get_loop_ranges - use element_mod, only: element_t - - type (element_t) , intent(inout), target :: elem(:) - integer , intent(in) :: nets,nete, nt - real(r8) :: dt - - real(r8) :: c1, c2 - integer :: k,ie - - do ie=nets,nete - do k=1,nlev - c2 = 1._r8 / (1._r8 + otau(k)*dt) - c1 = -otau(k) * c2 * dt - elem(ie)%state%v(:,:,:,k,nt) = elem(ie)%state%v(:,:,:,k,nt)+c1 * elem(ie)%state%v(:,:,:,k,nt) -! ptend%s(:ncol,k) = c3 * (state%u(:ncol,k)**2 + state%v(:ncol,k)**2) - enddo - end do - end subroutine rayleigh_friction - - - - subroutine solve_diffusion(dt,nx,nlev,i,j,nlay,pmid,pint,km,fld,boundary_condition,dfld) - use physconst, only: gravit - real(kind=r8), intent(in) :: dt - integer , intent(in) :: nlay, nlev,nx, i, j - real(kind=r8), intent(in) :: pmid(nx,nx,nlay),pint(nx,nx,nlay+1),km(nx,nx,nlay+1) - real(kind=r8), intent(in) :: fld(nx,nx,nlev) - real(kind=r8), intent(out) :: dfld(nlay) - integer :: boundary_condition - ! - real(kind=r8), dimension(nlay) :: current_guess,next_iterate - real(kind=r8) :: alp, alm, value_level0 - integer :: k,iter, niterations=4 - - ! Make the guess for the next time step equal to the initial value - current_guess(:)= fld(i,j,1:nlay) - do iter = 1, niterations - ! two formulations of the upper boundary condition - !next_iterate(1) = (initial_value(1) + alp * current_guess(i+1) + alm * current_guess(1)) /(1. + alp + alm) ! top BC, u'=0 - if (boundary_condition==0) then - next_iterate(1) = fld(i,j,1) ! u doesn't get prognosed by diffusion at top - else if (boundary_condition==1) then - value_level0 = 0.75_r8*fld(i,j,1) ! value above sponge - k=1 - alp = dt*(km(i,j,k+1)*gravit*gravit/(pmid(i,j,k)-pmid(i,j,k+1)))/(pint(i,j,k)-pint(i,j,k+1)) - alm = dt*(km(i,j,k )*gravit*gravit/(0.5_r8*(pmid(i,j,1)-pmid(i,j,2))))/(pint(i,j,k)-pint(i,j,k+1)) - next_iterate(k) = (fld(i,j,k) + alp * current_guess(k+1) + alm * value_level0)/(1._r8 + alp + alm) - else - ! - ! set fld'=0 at model top - ! - k=1 - alp = dt*(km(i,j,k+1)*gravit*gravit/(pmid(i,j,k)-pmid(i,j,k+1)))/(pint(i,j,k)-pint(i,j,k+1)) - alm = dt*(km(i,j,k )*gravit*gravit/(0.5_r8*(pmid(i,j,1)-pmid(i,j,2))))/(pint(i,j,k)-pint(i,j,k+1)) - next_iterate(k) = (fld(i,j,1) + alp * current_guess(2) + alm * current_guess(1))/(1._r8 + alp + alm) - end if - do k = 2, nlay-1 - alp = dt*(km(i,j,k+1)*gravit*gravit/(pmid(i,j,k )-pmid(i,j,k+1)))/(pint(i,j,k)-pint(i,j,k+1)) - alm = dt*(km(i,j,k )*gravit*gravit/(pmid(i,j,k-1)-pmid(i,j,k )))/(pint(i,j,k)-pint(i,j,k+1)) - next_iterate(k) = (fld(i,j,k) + alp * current_guess(k+1) + alm * current_guess(k-1))/(1._r8 + alp + alm) - end do - next_iterate(nlay) = (fld(i,j,nlay) + alp * fld(i,j,nlay) + alm * current_guess(nlay-1))/(1._r8 + alp + alm) ! bottom BC - - ! before the next iterate, make the current guess equal to the values of the last iteration - current_guess(:) = next_iterate(:) - end do - dfld(:) = next_iterate(:) - fld(i,j,1:nlay) - - end subroutine solve_diffusion - - end module prim_advance_mod diff --git a/src/dynamics/se/dycore/viscosity_mod.F90 b/src/dynamics/se/dycore/viscosity_mod.F90 index b29e48a1fa..714b731169 100644 --- a/src/dynamics/se/dycore/viscosity_mod.F90 +++ b/src/dynamics/se/dycore/viscosity_mod.F90 @@ -9,7 +9,7 @@ module viscosity_mod ! use shr_kind_mod, only: r8=>shr_kind_r8 use thread_mod, only: max_num_threads, omp_get_num_threads - use dimensions_mod, only: np, nc, nlev,qsize,nelemd + use dimensions_mod, only: np, nc, nlev,nlevp, qsize,nelemd use hybrid_mod, only: hybrid_t, get_loop_ranges, config_thread_region use parallel_mod, only: parallel_t use element_mod, only: element_t @@ -50,11 +50,11 @@ module viscosity_mod CONTAINS -subroutine biharmonic_wk_dp3d(elem,dptens,dpflux,ttens,vtens,deriv,edge3,hybrid,nt,nets,nete,kbeg,kend,& - dp3d_ref,T_ref) +subroutine biharmonic_wk_dp3d(elem,dptens,dpflux,ttens,vtens,deriv,edge3,hybrid,nt,nets,nete,kbeg,kend,hvcoord,& + dp3d_ref,T_ref,pmid_ref) use derivative_mod, only : subcell_Laplace_fluxes use dimensions_mod, only : ntrac, nu_div_lev,nu_lev - + use hybvcoord_mod, only : hvcoord_t !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! compute weak biharmonic operator ! input: h,v (stored in elem()%, in lat-lon coordinates @@ -68,17 +68,22 @@ subroutine biharmonic_wk_dp3d(elem,dptens,dpflux,ttens,vtens,deriv,edge3,hybrid, real (kind=r8), intent(out), dimension(nc,nc,4,nlev,nets:nete) :: dpflux real (kind=r8), dimension(np,np,2,nlev,nets:nete) :: vtens real (kind=r8), dimension(np,np,nlev,nets:nete) :: ttens,dptens - real (kind=r8), dimension(np,np,nlev,nets:nete), optional :: dp3d_ref,T_ref + real (kind=r8), dimension(np,np,nlev,nets:nete), optional :: dp3d_ref,T_ref,pmid_ref type (EdgeBuffer_t) , intent(inout) :: edge3 type (derivative_t) , intent(in) :: deriv - + type (hvcoord_t) , intent(in) :: hvcoord ! local integer :: i,j,k,kptr,ie,kblk ! real (kind=r8), dimension(:,:), pointer :: rspheremv real (kind=r8), dimension(np,np) :: tmp real (kind=r8), dimension(np,np) :: tmp2 real (kind=r8), dimension(np,np,2) :: v - real (kind=r8) :: nu_ratio1, nu_ratio2 + + real (kind=r8), dimension(np,np,nlev) :: lap_p_wk + real (kind=r8), dimension(np,np,nlevp) :: T_i + + + real (kind=r8) :: nu_ratio1, nu_ratio2, dp_thresh logical var_coef1 kblk = kend - kbeg + 1 @@ -88,8 +93,7 @@ subroutine biharmonic_wk_dp3d(elem,dptens,dpflux,ttens,vtens,deriv,edge3,hybrid, !so tensor is only used on second call to laplace_sphere_wk var_coef1 = .true. if(hypervis_scaling > 0) var_coef1 = .false. - - + dp_thresh=.025_r8 ! tunable coefficient do ie=nets,nete !$omp parallel do num_threads(vert_num_threads) private(k,tmp) do k=kbeg,kend @@ -123,7 +127,37 @@ subroutine biharmonic_wk_dp3d(elem,dptens,dpflux,ttens,vtens,deriv,edge3,hybrid, call vlaplace_sphere_wk(elem(ie)%state%v(:,:,:,k,nt),deriv,elem(ie),.true.,vtens(:,:,:,k,ie), & var_coef=var_coef1,nu_ratio=nu_ratio1) - enddo + enddo + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! add p correction to approximate Laplace on pressure surfaces + ! Laplace_p(T) = Laplace(T) - dT/dp Laplace(p) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + ! lap_p_wk should be precomputed: + do k=1,nlev + call laplace_sphere_wk(pmid_ref(:,:,k,ie),deriv,elem(ie),lap_p_wk(:,:,k),var_coef=.false.) + enddo + + ! average T to interfaces, then compute dT/dp on midpoints: + T_i(:,:,1) = elem(ie)%state%T(:,:,1,nt) + T_i(:,:,nlevp) = elem(ie)%state%T(:,:,nlev,nt) + do k=2,nlev + T_i(:,:,k)=(elem(ie)%state%T(:,:,k,nt) + elem(ie)%state%T(:,:,k-1,nt))/2 + enddo + + do k=1,nlev + if (hvcoord%hybm(k)>0) then + tmp(:,:) = (T_i(:,:,k+1)-T_i(:,:,k))/dp3d_ref(:,:,k,ie) + tmp(:,:)=tmp(:,:) / (1.0_r8 + abs(tmp(:,:))/dp_thresh) + ttens(:,:,k,ie)=ttens(:,:,k,ie)-tmp(:,:)*lap_p_wk(:,:,k) ! correction term + endif + enddo + + + + + kptr = kbeg - 1 call edgeVpack(edge3,ttens(:,:,kbeg:kend,ie),kblk,kptr,ie) diff --git a/src/dynamics/se/dyn_comp.F90 b/src/dynamics/se/dyn_comp.F90 index 231f445abd..e63be25ebc 100644 --- a/src/dynamics/se/dyn_comp.F90 +++ b/src/dynamics/se/dyn_comp.F90 @@ -115,7 +115,8 @@ subroutine dyn_readnl(NLFileName) use control_mod, only: topology, phys_dyn_cp, variable_nsplit use control_mod, only: fine_ne, hypervis_power, hypervis_scaling use control_mod, only: max_hypervis_courant, statediag_numtrac,refined_mesh - use control_mod, only: raytau0, raykrange, rayk0, molecular_diff + use control_mod, only: molecular_diff + use control_mod, only: sponge_del4_nu_div_fac, sponge_del4_nu_fac, sponge_del4_lev use dimensions_mod, only: ne, npart use dimensions_mod, only: lcp_moist use dimensions_mod, only: hypervis_dynamic_ref_state,large_Courant_incr @@ -153,6 +154,9 @@ subroutine dyn_readnl(NLFileName) real(r8) :: se_nu_div real(r8) :: se_nu_p real(r8) :: se_nu_top + real(r8) :: se_sponge_del4_nu_fac + real(r8) :: se_sponge_del4_nu_div_fac + integer :: se_sponge_del4_lev integer :: se_qsplit logical :: se_refined_mesh integer :: se_rsplit @@ -173,9 +177,6 @@ subroutine dyn_readnl(NLFileName) integer :: se_kmin_jet integer :: se_kmax_jet integer :: se_phys_dyn_cp - real(r8) :: se_raytau0 - real(r8) :: se_raykrange - integer :: se_rayk0 real(r8) :: se_molecular_diff namelist /dyn_se_inparm/ & @@ -198,6 +199,9 @@ subroutine dyn_readnl(NLFileName) se_nu_div, & se_nu_p, & se_nu_top, & + se_sponge_del4_nu_fac, & + se_sponge_del4_nu_div_fac, & + se_sponge_del4_lev, & se_qsplit, & se_refined_mesh, & se_rsplit, & @@ -221,9 +225,6 @@ subroutine dyn_readnl(NLFileName) se_kmin_jet, & se_kmax_jet, & se_phys_dyn_cp, & - se_raytau0, & - se_raykrange, & - se_rayk0, & se_molecular_diff !-------------------------------------------------------------------------- @@ -273,6 +274,9 @@ subroutine dyn_readnl(NLFileName) call MPI_bcast(se_nu_div, 1, mpi_real8, masterprocid, mpicom, ierr) call MPI_bcast(se_nu_p, 1, mpi_real8, masterprocid, mpicom, ierr) call MPI_bcast(se_nu_top, 1, mpi_real8, masterprocid, mpicom, ierr) + call MPI_bcast(se_sponge_del4_nu_fac, 1, mpi_real8, masterprocid, mpicom, ierr) + call MPI_bcast(se_sponge_del4_nu_div_fac, 1, mpi_real8, masterprocid, mpicom, ierr) + call MPI_bcast(se_sponge_del4_lev, 1, mpi_integer, masterprocid, mpicom, ierr) call MPI_bcast(se_qsplit, 1, mpi_integer, masterprocid, mpicom, ierr) call MPI_bcast(se_refined_mesh, 1, mpi_logical, masterprocid, mpicom, ierr) call MPI_bcast(se_rsplit, 1, mpi_integer, masterprocid, mpicom, ierr) @@ -297,9 +301,6 @@ subroutine dyn_readnl(NLFileName) call MPI_bcast(se_kmin_jet, 1, mpi_integer, masterprocid, mpicom, ierr) call MPI_bcast(se_kmax_jet, 1, mpi_integer, masterprocid, mpicom, ierr) call MPI_bcast(se_phys_dyn_cp, 1, mpi_integer, masterprocid, mpicom, ierr) - call MPI_bcast(se_rayk0 , 1, mpi_integer, masterprocid, mpicom, ierr) - call MPI_bcast(se_raykrange, 1, mpi_real8, masterprocid, mpicom, ierr) - call MPI_bcast(se_raytau0, 1, mpi_real8, masterprocid, mpicom, ierr) call MPI_bcast(se_molecular_diff, 1, mpi_real8, masterprocid, mpicom, ierr) if (se_npes <= 0) then @@ -351,6 +352,9 @@ subroutine dyn_readnl(NLFileName) nu_p = se_nu_p nu_q = se_nu_p !for tracer-wind consistency nu_q must me equal to nu_p nu_top = se_nu_top + sponge_del4_nu_fac = se_sponge_del4_nu_fac + sponge_del4_nu_div_fac = se_sponge_del4_nu_div_fac + sponge_del4_lev = se_sponge_del4_lev qsplit = se_qsplit rsplit = se_rsplit statefreq = se_statefreq @@ -367,9 +371,6 @@ subroutine dyn_readnl(NLFileName) kmax_jet = se_kmax_jet variable_nsplit = .false. phys_dyn_cp = se_phys_dyn_cp - raytau0 = se_raytau0 - raykrange = se_raykrange - rayk0 = se_rayk0 molecular_diff = se_molecular_diff if (fv_nphys > 0) then @@ -466,6 +467,14 @@ subroutine dyn_readnl(NLFileName) write(iulog, '(a,i0)') 'dyn_readnl: se_fvm_supercycling_jet = ',fvm_supercycling_jet write(iulog, '(a,i0)') 'dyn_readnl: se_kmin_jet = ',kmin_jet write(iulog, '(a,i0)') 'dyn_readnl: se_kmax_jet = ',kmax_jet + + write(iulog, *) 'dyn_readnl: se_sponge_del4_nu_fac = ',se_sponge_del4_nu_fac + if (se_sponge_del4_nu_fac < 0) write(iulog, '(a)') ' (automatically set based on model top location)' + write(iulog, *) 'dyn_readnl: se_sponge_del4_nu_div_fac = ',se_sponge_del4_nu_div_fac + if (se_sponge_del4_nu_div_fac < 0) write(iulog, '(a)') ' (automatically set based on model top location)' + write(iulog, *) 'dyn_readnl: se_sponge_del4_lev = ',se_sponge_del4_lev + if (se_sponge_del4_lev < 0) write(iulog, '(a)') ' (automatically set based on model top location)' + if (se_refined_mesh) then write(iulog, '(a)') 'dyn_readnl: Refined mesh simulation' write(iulog, '(a)') 'dyn_readnl: se_mesh_file = ',trim(se_mesh_file) @@ -499,9 +508,6 @@ subroutine dyn_readnl(NLFileName) write(iulog,'(a,l1)') 'dyn_readnl: write restart data on unstructured grid = ', & se_write_restart_unstruct - write(iulog, '(a,e9.2)') 'dyn_readnl: se_raytau0 = ', raytau0 - write(iulog, '(a,e9.2)') 'dyn_readnl: se_raykrange = ', raykrange - write(iulog, '(a,i0)' ) 'dyn_readnl: se_rayk0 = ', rayk0 write(iulog, '(a,e9.2)') 'dyn_readnl: se_molecular_diff = ', molecular_diff end if @@ -598,16 +604,16 @@ subroutine dyn_init(dyn_in, dyn_out) use hybrid_mod, only: get_loop_ranges, config_thread_region use dimensions_mod, only: nu_scale_top, nu_lev, nu_div_lev use dimensions_mod, only: ksponge_end, kmvis_ref, kmcnd_ref,rho_ref,km_sponge_factor - use dimensions_mod, only: kmvisi_ref, kmcndi_ref,rhoi_ref use dimensions_mod, only: cnst_name_gll, cnst_longname_gll - use dimensions_mod, only: irecons_tracer_lev, irecons_tracer, otau, kord_tr, kord_tr_cslam + use dimensions_mod, only: irecons_tracer_lev, irecons_tracer, kord_tr, kord_tr_cslam use prim_driver_mod, only: prim_init2 use time_mod, only: time_at - use control_mod, only: runtype, raytau0, raykrange, rayk0, molecular_diff, nu_top + use control_mod, only: runtype, molecular_diff, nu_top use test_fvm_mapping, only: test_mapping_addfld use phys_control, only: phys_getopts use physconst, only: get_molecular_diff_coef_reference use control_mod, only: vert_remap_uvTq_alg, vert_remap_tracer_alg + use std_atm_profile, only: std_atm_height ! Dummy arguments: type(dyn_import_t), intent(out) :: dyn_in type(dyn_export_t), intent(out) :: dyn_out @@ -615,7 +621,7 @@ subroutine dyn_init(dyn_in, dyn_out) ! Local variables integer :: ithr, nets, nete, ie, k, kmol_end real(r8), parameter :: Tinit = 300.0_r8 - real(r8) :: press, ptop, tref + real(r8) :: press(1), ptop, tref,z(1) type(hybrid_t) :: hybrid @@ -663,7 +669,6 @@ subroutine dyn_init(dyn_in, dyn_out) character(len=*), parameter :: subname = 'dyn_init' - real(r8) :: tau0, krange, otau0, scale real(r8) :: km_sponge_factor_local(nlev+1) !---------------------------------------------------------------------------- @@ -740,27 +745,6 @@ subroutine dyn_init(dyn_in, dyn_out) call clean_iodesc_list() end if ! - ! initialize Rayleigh friction - ! - krange = raykrange - if (raykrange .eq. 0._r8) krange = (rayk0 - 1) / 2._r8 - tau0 = (86400._r8) * raytau0 ! convert to seconds - otau0 = 0._r8 - if (tau0 .ne. 0._r8) otau0 = 1._r8/tau0 - do k = 1, nlev - otau(k) = otau0 * (1.0_r8 + tanh((rayk0 - k) / krange)) / (2._r8) - enddo - if (masterproc) then - if (tau0 > 0._r8) then - write (iulog,*) 'SE dycore Rayleigh friction - krange = ', krange - write (iulog,*) 'SE dycore Rayleigh friction - otau0 = ', 1.0_r8/tau0 - write (iulog,*) 'SE dycore Rayleigh friction decay rate profile (only applied to (u,v))' - do k = 1, nlev - write (iulog,*) ' k = ', k, ' otau = ', otau(k) - enddo - end if - end if - ! ! initialize diffusion in dycore ! kmol_end = 0 @@ -772,9 +756,6 @@ subroutine dyn_init(dyn_in, dyn_out) tref = 1000._r8 !mean value at model top for solar max km_sponge_factor = molecular_diff km_sponge_factor_local = molecular_diff - call get_molecular_diff_coef_reference(1,nlev+1,tref,& - (hvcoord%hyai(:)+hvcoord%hybi(:))*hvcoord%ps0, km_sponge_factor_local,& - kmvisi_ref,kmcndi_ref,rhoi_ref) ! ! get rho, kmvis and kmcnd at mid-levels ! @@ -786,8 +767,10 @@ subroutine dyn_init(dyn_in, dyn_out) ! only apply molecular viscosity where viscosity is > 1000 m/s^2 if (MIN(kmvis_ref(k)/rho_ref(k),kmcnd_ref(k)/(cpair*rho_ref(k)))>1000.0_r8) then if (masterproc) then - write(iulog,'(a,i3,2e11.4)') "k, p, km_sponge_factor :",k, & - (hvcoord%hyam(k)+hvcoord%hybm(k))*hvcoord%ps0,km_sponge_factor(k) + press = (hvcoord%hyam(k)+hvcoord%hybm(k))*hvcoord%ps0 + call std_atm_height(press,z) + write(iulog,'(a,i3,3e11.4)') "k, p, z, km_sponge_factor :",k, & + press, z,km_sponge_factor(k) write(iulog,'(a,2e11.4)') "kmvis_ref/rho_ref, kmcnd_ref/(cp*rho_ref): ", & kmvis_ref(k)/rho_ref(k),kmcnd_ref(k)/(cpair*rho_ref(k)) end if @@ -815,7 +798,7 @@ subroutine dyn_init(dyn_in, dyn_out) do k=1,nlev press = (hvcoord%hyam(k)+hvcoord%hybm(k))*hvcoord%ps0 ptop = hvcoord%hyai(1)*hvcoord%ps0 - nu_scale_top(k) = 8.0_r8*(1.0_r8+tanh(1.0_r8*log(ptop/press))) ! tau will be maximum 8 at model top + nu_scale_top(k) = 8.0_r8*(1.0_r8+tanh(1.0_r8*log(ptop/press(1)))) ! tau will be maximum 8 at model top if (nu_scale_top(k).ge.0.15_r8) then ksponge_end = k else @@ -829,8 +812,11 @@ subroutine dyn_init(dyn_in, dyn_out) if (masterproc) then write(iulog,*) subname//": ksponge_end = ",ksponge_end if (nu_top>0) then - do k=1,ksponge_end - write(iulog,'(a,i3,1e11.4)') subname//": nu_scale_top ",k,nu_scale_top(k) + do k=1,ksponge_end+1 + press = (hvcoord%hyam(k)+hvcoord%hybm(k))*hvcoord%ps0 + call std_atm_height(press,z) + write(iulog,'(a,i3,4e11.4)') subname//": k, p, z, nu_scale_top, nu ",k,press,z,& + nu_scale_top(k),nu_scale_top(k)*nu_top end do end if end if diff --git a/src/physics/cam/check_energy.F90 b/src/physics/cam/check_energy.F90 index 6fd157628f..4c71215aa9 100644 --- a/src/physics/cam/check_energy.F90 +++ b/src/physics/cam/check_energy.F90 @@ -985,6 +985,9 @@ subroutine calc_te_and_aam_budgets(state, outfld_name_suffix) mr_cnst = rearth**3/gravit mo_cnst = omega*rearth**4/gravit + + mr = 0.0_r8 + mo = 0.0_r8 do k = 1, pver do i = 1, ncol cos_lat = cos(state%lat(i)) @@ -996,7 +999,7 @@ subroutine calc_te_and_aam_budgets(state, outfld_name_suffix) end do end do call outfld(name_out1 ,mr, pcols,lchnk ) - call outfld(name_out1 ,mo, pcols,lchnk ) + call outfld(name_out2 ,mo, pcols,lchnk ) end if end subroutine calc_te_and_aam_budgets diff --git a/src/physics/cam/geopotential.F90 b/src/physics/cam/geopotential.F90 index 51c8b636b8..b06b145e51 100644 --- a/src/physics/cam/geopotential.F90 +++ b/src/physics/cam/geopotential.F90 @@ -77,7 +77,7 @@ subroutine geopotential_dse( & rog(:ncol,:) = rair(:ncol,:) / gravit ! set calculation method based on dycore type - calc1 = dycore_is ('LR').or.dycore_is ('SE').or. dycore_is('FV3') + calc1 = dycore_is ('LR').or.dycore_is('FV3') ! The surface height is zero by definition. do i = 1,ncol diff --git a/src/physics/cam/physpkg.F90 b/src/physics/cam/physpkg.F90 index 42818734de..0dd7c5c954 100644 --- a/src/physics/cam/physpkg.F90 +++ b/src/physics/cam/physpkg.F90 @@ -1860,7 +1860,8 @@ subroutine tphysac (ztodt, cam_in, & ! For not ('FV'|'FV3'), physics_dme_adjust is called for energy diagnostic purposes only. So, save off tracers if (.not.(dycore_is('FV').or.dycore_is('FV3')).and.& (hist_fld_active('SE_pAM').or.hist_fld_active('KE_pAM').or.hist_fld_active('WV_pAM').or.& - hist_fld_active('WL_pAM').or.hist_fld_active('WI_pAM'))) then + hist_fld_active('WL_pAM').or.hist_fld_active('WI_pAM').or.hist_fld_active('MR_pAM').or.& + hist_fld_active('MO_pAM'))) then tmp_trac(:ncol,:pver,:pcnst) = state%q(:ncol,:pver,:pcnst) tmp_pdel(:ncol,:pver) = state%pdel(:ncol,:pver) tmp_ps(:ncol) = state%ps(:ncol) diff --git a/src/physics/cam/tracers.F90 b/src/physics/cam/tracers.F90 index 3b1773745c..bd5dff976f 100644 --- a/src/physics/cam/tracers.F90 +++ b/src/physics/cam/tracers.F90 @@ -291,7 +291,11 @@ subroutine tracers_init_cnst(name, latvals, lonvals, mask, q, z) do m = 1, test_tracer_num if (name == test_tracer_names(m)) then if (analytic_tracer(m)) then - call test_func_set(name, latvals, lonvals, mask, q, z=z) + if (present(z)) then + call test_func_set(name, latvals, lonvals, mask, q, z=z) + else + call test_func_set(name, latvals, lonvals, mask, q) + end if found = .true. exit else @@ -528,10 +532,10 @@ function test_func(name, lat, lon, k, z) result(fout) fout = 2.0_r8 + cos(lon) case('TT_COSB') ! - ! Cosine bell (Kent et al., 2012, MWR) + ! Cosine bell inspired by Kent et al., 2012, MWR; only one bell and location changed ! https://journals.ametsoc.org/doi/pdf/10.1175/MWR-D-11-00150.1 ! - R0 = 0.9_r8*1.0_r8/2.0_r8 ! radius of the perturbation + R0 = 0.5_r8 ! radius of the perturbation lon1 = pi/9.0_r8 lat1 = 2.0_r8*pi/9.0_r8 @@ -546,18 +550,15 @@ function test_func(name, lat, lon, k, z) result(fout) end if if (Rg1 < R0) then - fout = 0.1_r8+0.5_r8*(1.0_r8+COS(pi*d1)) + fout = 0.1_r8+0.9_r8*0.5_r8*(1.0_r8+COS(pi*d1)) else fout = 0.1_r8 end if - ! IF (ABS(fout) < 1.0E-8_r8) fout = 0.0_r8 - ! eta_c = 0.6_r8 - ! eta = (hyam(k)*ps0 + hybm(k)*psurf_moist)/psurf_moist case('TT_CCOSB') ! ! Correlated cosine bell ! - R0 = 0.9_r8*1.0_r8/2.0_r8 ! radius of the perturbation + R0 = 0.5_r8 ! radius of the perturbation lon1 = pi/9.0_r8 lat1 = 2.0_r8*pi/9.0_r8 @@ -572,7 +573,7 @@ function test_func(name, lat, lon, k, z) result(fout) end if if (Rg1 < R0) then - f1 = 0.1_r8+0.5_r8*(1.0_r8+COS(pi*d1)) + f1 = 0.1_r8+0.9_r8*0.5_r8*(1.0_r8+COS(pi*d1)) else f1 = 0.1_r8 end if @@ -582,7 +583,7 @@ function test_func(name, lat, lon, k, z) result(fout) ! ! Correlated cosine bell ! - R0 = 0.9_r8*1.0_r8/2.0_r8 ! radius of the perturbation + R0 = 0.5_r8 ! radius of the perturbation lon1 = pi/9.0_r8 lat1 = 2.0_r8*pi/9.0_r8 @@ -597,7 +598,7 @@ function test_func(name, lat, lon, k, z) result(fout) end if if (Rg1 < R0) then - f1 = 0.1_r8+0.5_r8*(1.0_r8+COS(pi*d1)) + f1 = 0.1_r8+0.9_r8*0.5_r8*(1.0_r8+COS(pi*d1)) else f1 = 0.1_r8 end if diff --git a/src/physics/simple/physpkg.F90 b/src/physics/simple/physpkg.F90 index 0440b54e07..e0c21b5157 100644 --- a/src/physics/simple/physpkg.F90 +++ b/src/physics/simple/physpkg.F90 @@ -525,7 +525,7 @@ subroutine tphysac (ztodt, cam_in, cam_out, state, tend, pbuf) cldiceini = 0.0_r8 end if - call calc_te_and_aam_budgets(state, 'pAP') + !========================= ! Compute physics tendency @@ -536,6 +536,8 @@ subroutine tphysac (ztodt, cam_in, cam_out, state, tend, pbuf) call physics_update(state, ptend, ztodt, tend) end if + call calc_te_and_aam_budgets(state, 'pAP') + ! FV: convert dry-type mixing ratios to moist here because ! physics_dme_adjust assumes moist. This is done in p_d_coupling for ! other dynamics. Bundy, Feb 2004. @@ -563,8 +565,9 @@ subroutine tphysac (ztodt, cam_in, cam_out, state, tend, pbuf) ! For not 'FV'|'FV3', physics_dme_adjust is called for energy diagnostic purposes only. ! So, save off tracers if (.not.(dycore_is('FV').or.dycore_is('FV3')) .and. & - (hist_fld_active('SE_pAM').or.hist_fld_active('KE_pAM').or.hist_fld_active('WV_pAM').or.& - hist_fld_active('WL_pAM').or.hist_fld_active('WI_pAM'))) then + (hist_fld_active('SE_pAM').or.hist_fld_active('KE_pAM').or.hist_fld_active('WV_pAM').or.& + hist_fld_active('WL_pAM').or.hist_fld_active('WI_pAM').or.hist_fld_active('MR_pAM').or.& + hist_fld_active('MO_pAM'))) then tmp_trac(:ncol,:pver,:pcnst) = state%q(:ncol,:pver,:pcnst) tmp_pdel(:ncol,:pver) = state%pdel(:ncol,:pver) tmp_ps(:ncol) = state%ps(:ncol) diff --git a/src/utils/physconst.F90 b/src/utils/physconst.F90 index bc6dee9f82..1bcbceb6e0 100644 --- a/src/utils/physconst.F90 +++ b/src/utils/physconst.F90 @@ -699,7 +699,7 @@ subroutine physconst_update(mmr, t, lchnk, ncol, to_moist_factor) call get_molecular_diff_coef(1,ncol,1,1,pver,pver,t(:ncol,:),1,sponge_factor,kmvis(:ncol,:,lchnk), & kmcnd(:ncol,:,lchnk), pcnst, tracer=mmr(:ncol,:,:), fact=to_moist_fact(:ncol,:), & active_species_idx_dycore=thermodynamic_active_species_idx) - + cappav(:ncol,:,lchnk) = rairv(:ncol,:,lchnk)/cpairv(:ncol,:,lchnk) end subroutine physconst_update ! !**************************************************************************************************************** @@ -861,7 +861,7 @@ subroutine get_pmid_from_dp(i0,i1,j0,j1,k0,k1,dp,ptop,pmid,pint) pint_local(:,:,k) = dp(:,:,k-1)+pint_local(:,:,k-1) end do - if (dycore_is ('LR').or.dycore_is ('SE')) then + if (dycore_is ('LR').or.dycore_is ('FV3')) then do k=k0,k1 pmid(:,:,k) = dp(:,:,k)/(log(pint_local(:,:,k+1))-log(pint_local(:,:,k))) end do @@ -1016,7 +1016,7 @@ subroutine get_gz_given_dp_Tv_Rdry(i0,i1,j0,j1,nlev,dp,T_v,R_dry,phis,ptop,gz,pm ! integrate hydrostatic eqn ! gzh = phis - if (dycore_is ('LR').or.dycore_is ('SE')) then + if (dycore_is ('LR').or.dycore_is ('FV3')) then do k=nlev,1,-1 Rdry_tv(:,:) = R_dry(:,:,k)*T_v(:,:,k) gz(:,:,k) = gzh(:,:)+Rdry_tv(:,:)*(1.0_r8-pint(:,:,k)/pmid_local(:,:,k)) @@ -1816,6 +1816,32 @@ subroutine get_molecular_diff_coef(i0,i1,j0,j1,k1,nlev,temp,get_at_interfaces,sp end do end do else if (get_at_interfaces==0) then + do k=1,k1 + do j=j0,j1 + do i=i0,i1 + kmvis(i,j,k) = 0.0_r8 + kmcnd(i,j,k) = 0.0_r8 + residual = 1.0_r8 + do icnst=1,dry_air_species_num-1 + ispecies = idx_local(icnst) + mm = tracer(i,j,k,ispecies)*factor(i,j,k) + kmvis(i,j,k) = kmvis(i,j,k)+thermodynamic_active_species_kv(icnst)* & + thermodynamic_active_species_mwi(icnst)*mm + kmcnd(i,j,k) = kmcnd(i,j,k)+thermodynamic_active_species_kc(icnst)* & + thermodynamic_active_species_mwi(icnst)*mm + residual = residual - mm + end do + icnst=dry_air_species_num + kmvis(i,j,k) = kmvis(i,j,k)+thermodynamic_active_species_kv(icnst)* & + thermodynamic_active_species_mwi(icnst)*residual + kmcnd(i,j,k) = kmcnd(i,j,k)+thermodynamic_active_species_kc(icnst)* & + thermodynamic_active_species_mwi(icnst)*residual + + kmvis(i,j,k) = kmvis(i,j,k)*mbarv(i,j,k)*temp(i,j,k)**kv4*1.e-7_r8 + kmcnd(i,j,k) = kmcnd(i,j,k)*mbarv(i,j,k)*temp(i,j,k)**kc4*1.e-5_r8 + enddo + enddo + end do else call endrun('get_molecular_diff_coef: get_at_interfaces must be 0 or 1') end if