diff --git a/.gitmodules b/.gitmodules index f899639357..6c2dbe2e68 100644 --- a/.gitmodules +++ b/.gitmodules @@ -36,7 +36,7 @@ [submodule "atmos_phys"] path = src/atmos_phys url = https://github.com/ESCOMP/atmospheric_physics - fxtag = atmos_phys0_08_000 + fxtag = atmos_phys0_09_000 fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics diff --git a/bld/configure b/bld/configure index 138436315a..4577a906aa 100755 --- a/bld/configure +++ b/bld/configure @@ -2149,6 +2149,7 @@ sub write_filepath # in the list of filepaths. print $fh "$camsrcdir/src/physics/cam\n"; print $fh "$camsrcdir/src/atmos_phys/to_be_ccppized\n"; + print $fh "$camsrcdir/src/atmos_phys/phys_utils\n"; #Add the CCPP'ized subdirectories print $fh "$camsrcdir/src/atmos_phys/schemes/tropopause_find\n"; diff --git a/doc/ChangeLog b/doc/ChangeLog index e57931b4cd..84e989f6de 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,5 +1,80 @@ =============================================================== +Tag name: cam6_4_071 +Originator(s): mwaxmonsky +Date: 26 February 2025 +One-line Summary: PBL_utils atmospheric_physics integration +Github PR URL: https://github.com/ESCOMP/CAM/pull/1235 + +Purpose of changes (include the issue number and title text for each relevant GitHub issue): + +Closes #1235 - Integrates PBL now from atmospheric_physics to enable further +vertical diffusion CCPP-ization + +Describe any changes made to build system: N/A + +Describe any changes made to the namelist: N/A + +List any changes to the defaults for the boundary datasets: N/A + +Describe any substantial timing or memory changes: N/A + +Code reviewed by: nusbaume + +List all files eliminated: N/A + +List all files added and what they do: N/A + +List all existing files that have been modified, and describe the changes: + +M .gitmodules +M src/atmos_phys +- Updates atmos_phys submodule + +M bld/configure +- Add phys utils which now contains refactored PBL utils code + +M src/physics/cam/clubb_intr.F90 +M src/physics/cam/eddy_diff_cam.F90 +M src/physics/cam/hb_diff.F90 +M src/physics/cam/pbl_utils.F90 +M src/physics/cam/vertical_diffusion.F90 +- Updates old PBL references to new PBL interface API + +M src/physics/cam/physpkg.F90 +M src/physics/cam7/physpkg.F90 + - Removing `pbl_utils_init()` call as no longer needed + + +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. + +derecho/intel/aux_cam: + +ERP_Ln9.f09_f09_mg17.FCSD_HCO.derecho_intel.cam-outfrq9s (Overall: FAIL) +SMS_Ld1.f09_f09_mg17.FCHIST_GC.derecho_intel.cam-outfrq1d (Overall: DIFF) + - pre-existing failures due to HEMCO not having reproducible results (issues #1018 and #856) + +SMS_D_Ln9.f19_f19_mg17.FXHIST.derecho_intel.cam-outfrq9s_amie (Overall: FAIL) +SMS_D_Ln9_P1280x1.ne0CONUSne30x8_ne0CONUSne30x8_mt12.FCHIST.derecho_intel.cam-outfrq9s (Overall: FAIL) + - pre-existing failures due to build-namelist error requiring CLM/CTSM external update + +derecho/nvhpc/aux_cam: ALL PASS + +izumi/nag/aux_cam: ALL PASS + +izumi/gnu/aux_cam: ALL PASS + +CAM tag used for the baseline comparison tests if different than previous +tag: + +Summarize any changes to answers: b4b + +=============================================================== +=============================================================== + Tag name: cam6_4_070 Originator(s): patcal, nusbaume Date: 22 February 2025 diff --git a/src/atmos_phys b/src/atmos_phys index 89b628646b..7031edf10a 160000 --- a/src/atmos_phys +++ b/src/atmos_phys @@ -1 +1 @@ -Subproject commit 89b628646b1506f36b35e67038552f09fb0662e6 +Subproject commit 7031edf10a3b2b63b92573cf03277ce8d65c073e diff --git a/src/physics/cam/clubb_intr.F90 b/src/physics/cam/clubb_intr.F90 index 9bbed56277..ccd655d2ab 100644 --- a/src/physics/cam/clubb_intr.F90 +++ b/src/physics/cam/clubb_intr.F90 @@ -20,13 +20,14 @@ module clubb_intr use shr_kind_mod, only: r8=>shr_kind_r8 use ppgrid, only: pver, pverp, pcols, begchunk, endchunk use phys_control, only: phys_getopts - use physconst, only: cpair, gravit, rga, latvap, latice, zvir, rh2o, karman, pi + use physconst, only: cpair, gravit, rga, latvap, latice, zvir, rh2o, karman, pi, rair use air_composition, only: rairv, cpairv use cam_history_support, only: max_fieldname_len use spmd_utils, only: masterproc use constituents, only: pcnst, cnst_add - use pbl_utils, only: calc_ustar, calc_obklen + use atmos_phys_pbl_utils,only: calc_friction_velocity, calc_kinematic_heat_flux, calc_ideal_gas_rrho, & + calc_kinematic_water_vapor_flux, calc_kinematic_buoyancy_flux, calc_obukhov_length use ref_pres, only: top_lev => trop_cloud_top_lev #ifdef CLUBB_SGS @@ -3491,8 +3492,8 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & ubar = sqrt(state1%u(i,pver)**2+state1%v(i,pver)**2) if (ubar < 0.25_r8) ubar = 0.25_r8 - call calc_ustar( state1%t(i,pver), state1%pmid(i,pver), cam_in%wsx(i), cam_in%wsy(i), & - rrho(i), ustar ) + rrho(i) = calc_ideal_gas_rrho(rair, state1%t(i,pver), state1%pmid(i,pver)) + ustar = calc_friction_velocity(cam_in%wsx(i), cam_in%wsy(i), rrho(i)) upwp_sfc(i) = -state1%u(i,pver)*ustar**2/ubar vpwp_sfc(i) = -state1%v(i,pver)*ustar**2/ubar @@ -4706,13 +4707,13 @@ subroutine clubb_tend_cam( state, ptend_all, pbuf, hdtime, & enddo ! diagnose surface friction and obukhov length (inputs to diagnose PBL depth) - rrho(1:ncol) = (rga)*(state1%pdel(1:ncol,pver)/dz_g(1:ncol,pver)) - call calc_ustar( ncol, state1%t(1:ncol,pver), state1%pmid(1:ncol,pver), cam_in%wsx(1:ncol), cam_in%wsy(1:ncol), & - rrho(1:ncol), ustar2(1:ncol)) + rrho (1:ncol) = calc_ideal_gas_rrho(rair, state1%t(1:ncol,pver), state1%pmid(1:ncol,pver)) + ustar2 (1:ncol) = calc_friction_velocity(cam_in%wsx(1:ncol), cam_in%wsy(1:ncol), rrho(1:ncol)) ! use correct qflux from coupler - call calc_obklen( ncol, th(1:ncol,pver), thv(1:ncol,pver), cam_in%cflx(1:ncol,1), cam_in%shf(1:ncol), & - rrho(1:ncol), ustar2(1:ncol), kinheat(1:ncol), kinwat(1:ncol), kbfs(1:ncol), & - obklen(1:ncol)) + kinheat(1:ncol) = calc_kinematic_heat_flux(cam_in%shf(1:ncol), rrho(1:ncol), cpair) + kinwat (1:ncol) = calc_kinematic_water_vapor_flux(cam_in%cflx(1:ncol,1), rrho(1:ncol)) + kbfs (1:ncol) = calc_kinematic_buoyancy_flux(kinheat(1:ncol), zvir, th(1:ncol,pver), kinwat(1:ncol)) + obklen (1:ncol) = calc_obukhov_length(thv(1:ncol,pver), ustar2(1:ncol), gravit, karman, kbfs(1:ncol)) dummy2(:) = 0._r8 dummy3(:) = 0._r8 diff --git a/src/physics/cam/eddy_diff_cam.F90 b/src/physics/cam/eddy_diff_cam.F90 index 1742bf5038..05ff18d9cd 100644 --- a/src/physics/cam/eddy_diff_cam.F90 +++ b/src/physics/cam/eddy_diff_cam.F90 @@ -424,16 +424,16 @@ subroutine compute_eddy_diff( pbuf, lchnk , ! May. 2008. ! !-------------------------------------------------------------------- ! - use diffusion_solver, only: compute_vdiff - use cam_history, only: outfld - use phys_debug_util, only: phys_debug_col - use air_composition, only: cpairv - use pbl_utils, only: calc_ustar, austausch_atm - use error_messages, only: handle_errmsg - use coords_1d, only: Coords1D - use wv_saturation, only: qsat - use eddy_diff, only: trbintd, caleddy - use physics_buffer, only: pbuf_get_field + use diffusion_solver, only: compute_vdiff + use cam_history, only: outfld + use phys_debug_util, only: phys_debug_col + use air_composition, only: cpairv + use atmos_phys_pbl_utils, only: calc_eddy_flux_coefficient, calc_ideal_gas_rrho, calc_friction_velocity + use error_messages, only: handle_errmsg + use coords_1d, only: Coords1D + use wv_saturation, only: qsat + use eddy_diff, only: trbintd, caleddy + use physics_buffer, only: pbuf_get_field ! --------------- ! ! Input Variables ! @@ -670,10 +670,11 @@ subroutine compute_eddy_diff( pbuf, lchnk , ! I am using updated wind, here. ! Compute ustar - call calc_ustar( ncol, tfd(:ncol,pver), pmid(:ncol,pver), & - taux(:ncol) - ksrftms(:ncol) * ufd(:ncol,pver), & ! Zonal wind stress - tauy(:ncol) - ksrftms(:ncol) * vfd(:ncol,pver), & ! Meridional wind stress - rrho(:ncol), ustar(:ncol)) + rrho(:ncol) = calc_ideal_gas_rrho(rair, tfd(:ncol,pver), pmid(:ncol,pver)) + ustar(:ncol) = calc_friction_velocity(taux(:ncol) - ksrftms(:ncol) * ufd(:ncol,pver), & ! Zonal wind stress + tauy(:ncol) - ksrftms(:ncol) * vfd(:ncol,pver), & ! Meridional wind stress + rrho(:ncol)) + minpblh(:ncol) = 100.0_r8 * ustar(:ncol) ! By construction, 'minpblh' is larger than 1 [m] when 'ustar_min = 0.01'. ! Calculate (qt,sl,n2,s2,ri) from a given set of (t,qv,ql,qi,u,v) @@ -694,8 +695,12 @@ subroutine compute_eddy_diff( pbuf, lchnk , ! Get free atmosphere exchange coefficients. This 'kvf' is not used in UW moist PBL scheme if (use_kvf) then - call austausch_atm(pcols, ncol, pver, ntop_eddy, nbot_eddy, & - ml2, ri, s2, kvf ) + kvf(:ncol,:) = 0.0_r8 + do k = ntop_eddy, nbot_eddy-1 + do i = 1, ncol + kvf(i,k+1) = calc_eddy_flux_coefficient(ml2(k), ri(i, k), s2(i, k)) + end do + end do else kvf = 0._r8 end if diff --git a/src/physics/cam/hb_diff.F90 b/src/physics/cam/hb_diff.F90 index a3bb11a17d..81ad8ff7bf 100644 --- a/src/physics/cam/hb_diff.F90 +++ b/src/physics/cam/hb_diff.F90 @@ -11,7 +11,6 @@ module hb_diff ! Private methods: ! trbintd initializes time dependent variables ! pblintd initializes time dependent variables that depend pbl depth - ! austausch_atm computes free atmosphere exchange coefficients ! austausch_pbl computes pbl exchange coefficients ! !---------------------------Code history-------------------------------- @@ -150,13 +149,16 @@ subroutine compute_hb_diff(ncol , & ! !----------------------------------------------------------------------- - use pbl_utils, only: virtem, calc_ustar, calc_obklen, austausch_atm + use atmos_phys_pbl_utils, only: calc_virtual_temperature, calc_friction_velocity, calc_obukhov_length, & + calc_eddy_flux_coefficient, calc_ideal_gas_rrho, calc_kinematic_heat_flux, calc_kinematic_water_vapor_flux, & + calc_kinematic_buoyancy_flux + use physconst, only: zvir, rair, gravit, karman !------------------------------Arguments-------------------------------- ! ! Input arguments ! - integer, intent(in) :: ncol ! number of atmospheric columns + integer, intent(in) :: ncol ! number of atmospheric columns real(r8), intent(in) :: th(pcols,pver) ! potential temperature [K] real(r8), intent(in) :: t(pcols,pver) ! temperature (used for density) @@ -203,20 +205,22 @@ subroutine compute_hb_diff(ncol , & real(r8) :: n2(pcols,pver) ! brunt vaisaila frequency real(r8) :: bge(pcols) ! buoyancy gradient enhancment integer :: ktopbl(pcols) ! index of first midpoint inside pbl + integer :: i,k ! ! Initialize time dependent variables that do not depend on pbl height ! ! virtual temperature - call virtem(ncol, (pver-ntop_turb+1), th(:ncol,ntop_turb:),q(:ncol,ntop_turb:), thv(:ncol,ntop_turb:)) + thv(:ncol,ntop_turb:) = calc_virtual_temperature(th(:ncol,ntop_turb:),q(:ncol,ntop_turb:), zvir) ! Compute ustar, Obukhov length, and kinematic surface fluxes. - call calc_ustar(ncol, t(:ncol,pver),pmid(:ncol,pver),taux(:ncol),tauy(:ncol), & - rrho(:ncol),ustar(:ncol)) - call calc_obklen(ncol, th(:ncol,pver), thv(:ncol,pver), qflx(:ncol), & - shflx(:ncol), rrho(:ncol), ustar(:ncol), & - khfs(:ncol), kqfs(:ncol), kbfs(:ncol), & - obklen(:ncol)) + rrho(:ncol) = calc_ideal_gas_rrho(rair, t(:ncol,pver), pmid(:ncol,pver)) + ustar(:ncol) = calc_friction_velocity(taux(:ncol),tauy(:ncol), rrho(:ncol)) + khfs(:ncol) = calc_kinematic_heat_flux(shflx(:ncol), rrho(:ncol), cpair) + kqfs(:ncol) = calc_kinematic_water_vapor_flux(qflx(:ncol), rrho(:ncol)) + kbfs(:ncol) = calc_kinematic_buoyancy_flux(khfs(:ncol), zvir, th(:ncol,pver), kqfs(:ncol)) + obklen(:ncol) = calc_obukhov_length(thv(:ncol,pver), ustar(:ncol), gravit, karman, kbfs(:ncol)) + ! Calculate s2, n2, and Richardson number. call trbintd(ncol , & thv ,z ,u ,v , & @@ -229,10 +233,15 @@ subroutine compute_hb_diff(ncol , & ustar ,obklen ,kbfs ,pblh ,wstar , & zi ,cldn ,ocnfrac ,bge ) ! - ! Get free atmosphere exchange coefficients + ! Get atmosphere exchange coefficients ! - call austausch_atm(pcols, ncol, pver, ntop_turb, nbot_turb, & - ml2, ri, s2, kvf) + kvf(:ncol,:) = 0.0_r8 + do k = ntop_turb, nbot_turb-1 + do i = 1, ncol + kvf(i,k+1) = calc_eddy_flux_coefficient(ml2(k), ri(i, k), s2(i, k)) + end do + end do + ! ! Get pbl exchange coefficients ! @@ -262,7 +271,10 @@ subroutine compute_hb_free_atm_diff(ncol, & ! !----------------------------------------------------------------------- - use pbl_utils, only: virtem, calc_ustar, calc_obklen, austausch_atm_free + use atmos_phys_pbl_utils, only: calc_virtual_temperature, calc_friction_velocity, calc_obukhov_length, & + calc_free_atm_eddy_flux_coefficient, calc_ideal_gas_rrho, calc_kinematic_heat_flux, calc_kinematic_water_vapor_flux, & + calc_kinematic_buoyancy_flux + use physconst, only: zvir, rair, gravit, karman !------------------------------Arguments-------------------------------- ! @@ -303,17 +315,18 @@ subroutine compute_hb_free_atm_diff(ncol, & real(r8) :: kvf(pcols,pverp) ! free atmospheric eddy diffsvty [m2/s] real(r8) :: s2(pcols,pver) ! shear squared real(r8) :: n2(pcols,pver) ! brunt vaisaila frequency + integer :: i, k ! virtual potential temperature - call virtem(ncol, (pver-ntop_turb+1), th(:ncol,ntop_turb:),q(:ncol,ntop_turb:), thv(:ncol,ntop_turb:)) + thv(:ncol,ntop_turb:) = calc_virtual_temperature(th(:ncol,ntop_turb:),q(:ncol,ntop_turb:), zvir) ! Compute ustar, Obukhov length, and kinematic surface fluxes. - call calc_ustar(ncol, t(:ncol,pver),pmid(:ncol,pver),taux(:ncol),tauy(:ncol), & - rrho(:ncol),ustar(:ncol)) - call calc_obklen(ncol, th(:ncol,pver), thv(:ncol,pver), qflx(:ncol), & - shflx(:ncol), rrho(:ncol), ustar(:ncol), & - khfs(:ncol), kqfs(:ncol), kbfs(:ncol), & - obklen(:ncol)) + rrho(:ncol) = calc_ideal_gas_rrho(rair, t(:ncol,pver), pmid(:ncol,pver)) + ustar(:ncol) = calc_friction_velocity(taux(:ncol),tauy(:ncol), rrho(:ncol)) + khfs(:ncol) = calc_kinematic_heat_flux(shflx(:ncol), rrho(:ncol), cpair) + kqfs(:ncol) = calc_kinematic_water_vapor_flux(qflx(:ncol), rrho(:ncol)) + kbfs(:ncol) = calc_kinematic_buoyancy_flux(khfs(:ncol), zvir, th(:ncol,pver), kqfs(:ncol)) + obklen(:ncol) = calc_obukhov_length(thv(:ncol,pver), ustar(:ncol), gravit, karman, kbfs(:ncol)) ! Calculate s2, n2, and Richardson number. call trbintd(ncol , & thv ,z ,u ,v , & @@ -321,8 +334,12 @@ subroutine compute_hb_free_atm_diff(ncol, & ! ! Get free atmosphere exchange coefficients ! - call austausch_atm_free(pcols, ncol, pver, ntop_turb, nbot_turb, & - ml2, ri, s2, kvf) + kvf(:ncol,:) = 0.0_r8 + do k = ntop_turb, nbot_turb - 1 + do i = 1, ncol + kvf(i,k+1) = calc_free_atm_eddy_flux_coefficient(ml2(k), ri(i, k), s2(i, k)) + end do + end do kvq(:ncol,:) = kvf(:ncol,:) kvm(:ncol,:) = kvf(:ncol,:) diff --git a/src/physics/cam/pbl_utils.F90 b/src/physics/cam/pbl_utils.F90 index 66759e295d..cd3af2ec60 100644 --- a/src/physics/cam/pbl_utils.F90 +++ b/src/physics/cam/pbl_utils.F90 @@ -22,206 +22,10 @@ module pbl_utils ! procedures, so they can accept scalars or any dimension of array as ! arguments, as long as all arguments have the same number of ! elements. -public pbl_utils_init -public calc_ustar -public calc_obklen -public virtem public compute_radf -public austausch_atm, austausch_atm_free - -real(r8), parameter :: ustar_min = 0.01_r8 - -real(r8) :: g ! acceleration of gravity -real(r8) :: vk ! Von Karman's constant -real(r8) :: cpair ! specific heat of dry air -real(r8) :: rair ! gas constant for dry air -real(r8) :: zvir ! rh2o/rair - 1 - - -!------------------------------------------------------------------------! -! Purpose: Compilers aren't creating optimized vector versions of ! -! elemental routines, so we'll explicitly create them and bind ! -! them via an interface for transparent use ! -!------------------------------------------------------------------------! -interface calc_ustar - module procedure calc_ustar_scalar - module procedure calc_ustar_vector -end interface - -interface calc_obklen - module procedure calc_obklen_scalar - module procedure calc_obklen_vector -end interface - -interface virtem - module procedure virtem_vector1D - module procedure virtem_vector2D ! Used in hb_diff.F90 -end interface - - contains -subroutine pbl_utils_init(g_in,vk_in,cpair_in,rair_in,zvir_in) - - !-----------------------------------------------------------------------! - ! Purpose: Set constants to be used in calls to later functions ! - !-----------------------------------------------------------------------! - - real(r8), intent(in) :: g_in ! acceleration of gravity - real(r8), intent(in) :: vk_in ! Von Karman's constant - real(r8), intent(in) :: cpair_in ! specific heat of dry air - real(r8), intent(in) :: rair_in ! gas constant for dry air - real(r8), intent(in) :: zvir_in ! rh2o/rair - 1 - - g = g_in - vk = vk_in - cpair = cpair_in - rair = rair_in - zvir = zvir_in - -end subroutine pbl_utils_init - -subroutine calc_ustar_scalar( t, pmid, taux, tauy, & - rrho, ustar) - - !-----------------------------------------------------------------------! - ! Purpose: Calculate ustar and bottom level density (necessary for ! - ! Obukhov length calculation). ! - !-----------------------------------------------------------------------! - - real(r8), intent(in) :: t ! surface temperature - real(r8), intent(in) :: pmid ! midpoint pressure (bottom level) - real(r8), intent(in) :: taux ! surface u stress [N/m2] - real(r8), intent(in) :: tauy ! surface v stress [N/m2] - - real(r8), intent(out) :: rrho ! 1./bottom level density - real(r8), intent(out) :: ustar ! surface friction velocity [m/s] - - rrho = rair * t / pmid - ustar = max( sqrt( sqrt(taux**2 + tauy**2)*rrho ), ustar_min ) - -end subroutine calc_ustar_scalar - -subroutine calc_ustar_vector(n, t, pmid, taux, tauy, & - rrho, ustar) - - !-----------------------------------------------------------------------! - ! Purpose: Calculate ustar and bottom level density (necessary for ! - ! Obukhov length calculation). ! - !-----------------------------------------------------------------------! - integer, intent(in) :: n ! Length of vectors - - real(r8), intent(in) :: t(n) ! surface temperature - real(r8), intent(in) :: pmid(n) ! midpoint pressure (bottom level) - real(r8), intent(in) :: taux(n) ! surface u stress [N/m2] - real(r8), intent(in) :: tauy(n) ! surface v stress [N/m2] - - - real(r8), intent(out) :: rrho(n) ! 1./bottom level density - real(r8), intent(out) :: ustar(n) ! surface friction velocity [m/s] - - - rrho = rair * t / pmid - ustar = max( sqrt( sqrt(taux**2 + tauy**2)*rrho ), ustar_min ) - -end subroutine calc_ustar_vector - -subroutine calc_obklen_scalar( ths, thvs, qflx, shflx, rrho, ustar, & - khfs, kqfs, kbfs, obklen) - - !-----------------------------------------------------------------------! - ! Purpose: Calculate Obukhov length and kinematic fluxes. ! - !-----------------------------------------------------------------------! - - real(r8), intent(in) :: ths ! potential temperature at surface [K] - real(r8), intent(in) :: thvs ! virtual potential temperature at surface - real(r8), intent(in) :: qflx ! water vapor flux (kg/m2/s) - real(r8), intent(in) :: shflx ! surface heat flux (W/m2) - - real(r8), intent(in) :: rrho ! 1./bottom level density [ m3/kg ] - real(r8), intent(in) :: ustar ! Surface friction velocity [ m/s ] - - real(r8), intent(out) :: khfs ! sfc kinematic heat flux [mK/s] - real(r8), intent(out) :: kqfs ! sfc kinematic water vapor flux [m/s] - real(r8), intent(out) :: kbfs ! sfc kinematic buoyancy flux [m^2/s^3] - real(r8), intent(out) :: obklen ! Obukhov length - - ! Need kinematic fluxes for Obukhov: - khfs = shflx*rrho/cpair - kqfs = qflx*rrho - kbfs = khfs + zvir*ths*kqfs - - ! Compute Obukhov length: - obklen = -thvs * ustar**3 / (g*vk*(kbfs + sign(1.e-10_r8,kbfs))) - -end subroutine calc_obklen_scalar - -subroutine calc_obklen_vector(n, ths, thvs, qflx, shflx, rrho, ustar, & - khfs, kqfs, kbfs, obklen) - - !-----------------------------------------------------------------------! - ! Purpose: Calculate Obukhov length and kinematic fluxes. ! - !-----------------------------------------------------------------------! - integer, intent(in) :: n ! Length of vectors - - real(r8), intent(in) :: ths(n) ! potential temperature at surface [K] - real(r8), intent(in) :: thvs(n) ! virtual potential temperature at surface - real(r8), intent(in) :: qflx(n) ! water vapor flux (kg/m2/s) - real(r8), intent(in) :: shflx(n) ! surface heat flux (W/m2) - - real(r8), intent(in) :: rrho(n) ! 1./bottom level density [ m3/kg ] - real(r8), intent(in) :: ustar(n) ! Surface friction velocity [ m/s ] - - real(r8), intent(out) :: khfs(n) ! sfc kinematic heat flux [mK/s] - real(r8), intent(out) :: kqfs(n) ! sfc kinematic water vapor flux [m/s] - real(r8), intent(out) :: kbfs(n) ! sfc kinematic buoyancy flux [m^2/s^3] - real(r8), intent(out) :: obklen(n) ! Obukhov length - - - ! Need kinematic fluxes for Obukhov: - khfs = shflx*rrho/cpair - kqfs = qflx*rrho - kbfs = khfs + zvir*ths*kqfs - - ! Compute Obukhov length: - obklen = -thvs * ustar**3 / (g*vk*(kbfs + sign(1.e-10_r8,kbfs))) - -end subroutine calc_obklen_vector - -subroutine virtem_vector1D(n, t,q, virtem) - - !-----------------------------------------------------------------------! - ! Purpose: Calculate virtual temperature from temperature and specific ! - ! humidity. ! - !-----------------------------------------------------------------------! - - integer, intent(in) :: n ! vector length - - real(r8), intent(in) :: t(n), q(n) - real(r8), intent(out):: virtem(n) - - virtem = t * (1.0_r8 + zvir*q) - -end subroutine virtem_vector1D - -subroutine virtem_vector2D(n, m, t, q, virtem) - - !-----------------------------------------------------------------------! - ! Purpose: Calculate virtual temperature from temperature and specific ! - ! humidity. ! - !-----------------------------------------------------------------------! - - integer, intent(in) :: n, m ! vector lengths - - real(r8), intent(in) :: t(n,m), q(n,m) - real(r8), intent(out):: virtem(n,m) - - virtem = t * (1.0_r8 + zvir*q) - -end subroutine virtem_vector2D - - subroutine compute_radf( choice_radf, i, pcols, pver, ncvmax, ncvfin, ktop, qmin, & ql, pi, qrlw, g, cldeff, zi, chs, lwp_CL, opt_depth_CL, & radinvfrac_CL, radf_CL ) @@ -332,138 +136,4 @@ subroutine compute_radf( choice_radf, i, pcols, pver, ncvmax, ncvfin, ktop, qmin end do ! ncv = 1, ncvfin(i) end subroutine compute_radf -subroutine austausch_atm(pcols, ncol, pver, ntop, nbot, ml2, ri, s2, kvf) - - !---------------------------------------------------------------------- ! - ! ! - ! Purpose: Computes exchange coefficients for free turbulent flows. ! - ! ! - ! Method: ! - ! ! - ! The free atmosphere diffusivities are based on standard mixing length ! - ! forms for the neutral diffusivity multiplied by functns of Richardson ! - ! number. K = l^2 * |dV/dz| * f(Ri). The same functions are used for ! - ! momentum, potential temperature, and constitutents. ! - ! ! - ! The stable Richardson num function (Ri>0) is taken from Holtslag and ! - ! Beljaars (1989), ECMWF proceedings. f = 1 / (1 + 10*Ri*(1 + 8*Ri)) ! - ! The unstable Richardson number function (Ri<0) is taken from CCM1. ! - ! f = sqrt(1 - 18*Ri) ! - ! ! - ! Author: B. Stevens (rewrite, August 2000) ! - ! ! - !---------------------------------------------------------------------- ! - - ! --------------- ! - ! Input arguments ! - ! --------------- ! - - integer, intent(in) :: pcols ! Atmospheric columns dimension size - integer, intent(in) :: ncol ! Number of atmospheric columns - integer, intent(in) :: pver ! Number of atmospheric layers - integer, intent(in) :: ntop ! Top layer for calculation - integer, intent(in) :: nbot ! Bottom layer for calculation - - real(r8), intent(in) :: ml2(pver+1) ! Mixing lengths squared - real(r8), intent(in) :: s2(pcols,pver) ! Shear squared - real(r8), intent(in) :: ri(pcols,pver) ! Richardson no - - ! ---------------- ! - ! Output arguments ! - ! ---------------- ! - - real(r8), intent(out) :: kvf(pcols,pver+1) ! Eddy diffusivity for heat and tracers - - ! --------------- ! - ! Local Variables ! - ! --------------- ! - - real(r8) :: fofri ! f(ri) - real(r8) :: kvn ! Neutral Kv - - integer :: i ! Longitude index - integer :: k ! Vertical index - - real(r8), parameter :: zkmin = 0.01_r8 ! Minimum kneutral*f(ri). - - ! ----------------------- ! - ! Main Computation Begins ! - ! ----------------------- ! - - kvf(:ncol,:) = 0.0_r8 - - ! Compute the free atmosphere vertical diffusion coefficients: kvh = kvq = kvm. - - do k = ntop, nbot - 1 - do i = 1, ncol - if( ri(i,k) < 0.0_r8 ) then - fofri = sqrt( max( 1._r8 - 18._r8 * ri(i,k), 0._r8 ) ) - else - fofri = 1.0_r8 / ( 1.0_r8 + 10.0_r8 * ri(i,k) * ( 1.0_r8 + 8.0_r8 * ri(i,k) ) ) - end if - kvn = ml2(k) * sqrt(s2(i,k)) - kvf(i,k+1) = max( zkmin, kvn * fofri ) - end do - end do - -end subroutine austausch_atm - -subroutine austausch_atm_free(pcols, ncol, pver, ntop, nbot, ml2, ri, s2, kvf) - - !---------------------------------------------------------------------- ! - ! ! - ! same as austausch_atm but only mixing for Ri<0 ! - ! i.e. no background mixing and mixing for Ri>0 ! - ! ! - !---------------------------------------------------------------------- ! - - ! --------------- ! - ! Input arguments ! - ! --------------- ! - - integer, intent(in) :: pcols ! Atmospheric columns dimension size - integer, intent(in) :: ncol ! Number of atmospheric columns - integer, intent(in) :: pver ! Number of atmospheric layers - integer, intent(in) :: ntop ! Top layer for calculation - integer, intent(in) :: nbot ! Bottom layer for calculation - - real(r8), intent(in) :: ml2(pver+1) ! Mixing lengths squared - real(r8), intent(in) :: s2(pcols,pver) ! Shear squared - real(r8), intent(in) :: ri(pcols,pver) ! Richardson no - - ! ---------------- ! - ! Output arguments ! - ! ---------------- ! - - real(r8), intent(out) :: kvf(pcols,pver+1) ! Eddy diffusivity for heat and tracers - - ! --------------- ! - ! Local Variables ! - ! --------------- ! - - real(r8) :: fofri ! f(ri) - real(r8) :: kvn ! Neutral Kv - - integer :: i ! Longitude index - integer :: k ! Vertical index - - ! ----------------------- ! - ! Main Computation Begins ! - ! ----------------------- ! - - kvf(:ncol,:) = 0.0_r8 - ! Compute the free atmosphere vertical diffusion coefficients: kvh = kvq = kvm. - do k = ntop, nbot - 1 - do i = 1, ncol - if( ri(i,k) < 0.0_r8 ) then - fofri = sqrt( max( 1._r8 - 18._r8 * ri(i,k), 0._r8 ) ) - else - fofri = 0.0_r8 - end if - kvn = ml2(k) * sqrt(s2(i,k)) - kvf(i,k+1) = kvn * fofri - end do - end do -end subroutine austausch_atm_free - end module pbl_utils diff --git a/src/physics/cam/physpkg.F90 b/src/physics/cam/physpkg.F90 index 0e83ad2707..577b45df58 100644 --- a/src/physics/cam/physpkg.F90 +++ b/src/physics/cam/physpkg.F90 @@ -742,7 +742,6 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) use tracers, only: tracers_init use aoa_tracers, only: aoa_tracers_init use rayleigh_friction, only: rayleigh_friction_init - use pbl_utils, only: pbl_utils_init use vertical_diffusion, only: vertical_diffusion_init use phys_debug_util, only: phys_debug_init use rad_constituents, only: rad_cnst_init @@ -877,7 +876,6 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) call rayleigh_friction_init() - call pbl_utils_init(gravit, karman, cpair, rair, zvir) call vertical_diffusion_init(pbuf2d) if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then diff --git a/src/physics/cam/vertical_diffusion.F90 b/src/physics/cam/vertical_diffusion.F90 index efd571269a..d44f952559 100644 --- a/src/physics/cam/vertical_diffusion.F90 +++ b/src/physics/cam/vertical_diffusion.F90 @@ -710,28 +710,30 @@ subroutine vertical_diffusion_tend( & use physics_buffer, only : physics_buffer_desc, pbuf_get_field, pbuf_set_field use physics_types, only : physics_state, physics_ptend, physics_ptend_init - use camsrfexch, only : cam_in_t - use cam_history, only : outfld - - use trb_mtn_stress_cam, only : trb_mtn_stress_tend - use beljaars_drag_cam, only : beljaars_drag_tend - use eddy_diff_cam, only : eddy_diff_tend - use hb_diff, only : compute_hb_diff, compute_hb_free_atm_diff - use wv_saturation, only : qsat - use molec_diff, only : compute_molec_diff, vd_lu_qdecomp - use constituents, only : qmincg, qmin, cnst_type - use diffusion_solver, only : compute_vdiff, any, operator(.not.) - use air_composition, only : cpairv, rairv !Needed for calculation of upward H flux - use time_manager, only : get_nstep - use constituents, only : cnst_get_type_byind, cnst_name, & - cnst_mw, cnst_fixed_ubc, cnst_fixed_ubflx - use physconst, only : pi - use pbl_utils, only : virtem, calc_obklen, calc_ustar - use upper_bc, only : ubc_get_vals, ubc_fixed_temp - use upper_bc, only : ubc_get_flxs - use coords_1d, only : Coords1D - use phys_control, only : cam_physpkg_is - use ref_pres, only : ptop_ref + use camsrfexch, only : cam_in_t + use cam_history, only : outfld + + use trb_mtn_stress_cam, only : trb_mtn_stress_tend + use beljaars_drag_cam, only : beljaars_drag_tend + use eddy_diff_cam, only : eddy_diff_tend + use hb_diff, only : compute_hb_diff, compute_hb_free_atm_diff + use wv_saturation, only : qsat + use molec_diff, only : compute_molec_diff, vd_lu_qdecomp + use constituents, only : qmincg, qmin, cnst_type + use diffusion_solver, only : compute_vdiff, any, operator(.not.) + use air_composition, only : cpairv, rairv !Needed for calculation of upward H flux + use time_manager, only : get_nstep + use constituents, only : cnst_get_type_byind, cnst_name, & + cnst_mw, cnst_fixed_ubc, cnst_fixed_ubflx + use physconst, only : pi + use atmos_phys_pbl_utils, only: calc_virtual_temperature, calc_ideal_gas_rrho, calc_friction_velocity, & + calc_kinematic_heat_flux, calc_kinematic_water_vapor_flux, calc_kinematic_buoyancy_flux, & + calc_obukhov_length + use upper_bc, only : ubc_get_vals, ubc_fixed_temp + use upper_bc, only : ubc_get_flxs + use coords_1d, only : Coords1D + use phys_control, only : cam_physpkg_is + use ref_pres, only : ptop_ref ! --------------- ! ! Input Arguments ! @@ -1020,11 +1022,12 @@ subroutine vertical_diffusion_tend( & ! The diag_TKE scheme does not calculate the Monin-Obukhov length, which is used in dry deposition calculations. ! Use the routines from pbl_utils to accomplish this. Assumes ustar and rrho have been set. - call virtem(ncol, th(:ncol,pver),state%q(:ncol,pver,1), thvs(:ncol)) - call calc_obklen(ncol, th(:ncol,pver), thvs(:ncol), cam_in%cflx(:ncol,1), & - cam_in%shf(:ncol), rrho(:ncol), ustar(:ncol), & - khfs(:ncol), kqfs(:ncol), kbfs(:ncol), obklen(:ncol)) + thvs (:ncol) = calc_virtual_temperature(th(:ncol,pver), state%q(:ncol,pver,1), zvir) + khfs (:ncol) = calc_kinematic_heat_flux(cam_in%shf(:ncol), rrho(:ncol), cpair) + kqfs (:ncol) = calc_kinematic_water_vapor_flux(cam_in%cflx(:ncol,1), rrho(:ncol)) + kbfs (:ncol) = calc_kinematic_buoyancy_flux(khfs(:ncol), zvir, th(:ncol,pver), kqfs(:ncol)) + obklen(:ncol) = calc_obukhov_length(thvs(:ncol), ustar(:ncol), gravit, karman, kbfs(:ncol)) case ( 'HB', 'HBR' ) @@ -1076,14 +1079,14 @@ subroutine vertical_diffusion_tend( & ! is only handling other things, e.g. some boundary conditions, tms, ! and molecular diffusion. - call virtem(ncol, th(:ncol,pver),state%q(:ncol,pver,1), thvs(:ncol)) + thvs (:ncol) = calc_virtual_temperature(th(:ncol,pver), state%q(:ncol,pver,1), zvir) + rrho (:ncol) = calc_ideal_gas_rrho(rair, state%t(:ncol,pver), state%pmid(:ncol,pver)) + ustar (:ncol) = calc_friction_velocity(cam_in%wsx(:ncol), cam_in%wsy(:ncol), rrho(:ncol)) + khfs (:ncol) = calc_kinematic_heat_flux(cam_in%shf(:ncol), rrho(:ncol), cpair) + kqfs (:ncol) = calc_kinematic_water_vapor_flux(cam_in%cflx(:ncol,1), rrho(:ncol)) + kbfs (:ncol) = calc_kinematic_buoyancy_flux(khfs(:ncol), zvir, th(:ncol,pver), kqfs(:ncol)) + obklen(:ncol) = calc_obukhov_length(thvs(:ncol), ustar(:ncol), gravit, karman, kbfs(:ncol)) - call calc_ustar( ncol, state%t(:ncol,pver), state%pmid(:ncol,pver), & - cam_in%wsx(:ncol), cam_in%wsy(:ncol), rrho(:ncol), ustar(:ncol)) - ! Use actual qflux, not lhf/latvap as was done previously - call calc_obklen( ncol, th(:ncol,pver), thvs(:ncol), cam_in%cflx(:ncol,1), & - cam_in%shf(:ncol), rrho(:ncol), ustar(:ncol), & - khfs(:ncol), kqfs(:ncol), kbfs(:ncol), obklen(:ncol)) ! These tendencies all applied elsewhere. kvm = 0._r8 kvh = 0._r8 diff --git a/src/physics/cam7/physpkg.F90 b/src/physics/cam7/physpkg.F90 index 43b55137f9..ffcd05293e 100644 --- a/src/physics/cam7/physpkg.F90 +++ b/src/physics/cam7/physpkg.F90 @@ -744,7 +744,6 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) use tracers, only: tracers_init use aoa_tracers, only: aoa_tracers_init use rayleigh_friction, only: rayleigh_friction_init - use pbl_utils, only: pbl_utils_init use vertical_diffusion, only: vertical_diffusion_init use phys_debug_util, only: phys_debug_init use phys_debug, only: phys_debug_state_init @@ -880,7 +879,6 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) call rayleigh_friction_init() - call pbl_utils_init(gravit, karman, cpair, rair, zvir) call vertical_diffusion_init(pbuf2d) if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then