From 1081b0a6ed60ea072a8d8230299c8d588f6574ea Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Wed, 13 May 2026 14:45:15 -0400 Subject: [PATCH 1/3] Initial CAM work (surprisingly hand-written) --- src/physics/cam/beljaars_drag.F90 | 235 +++++++++++--------------- src/physics/cam/beljaars_drag_cam.F90 | 48 ++++-- src/physics/cam/physpkg.F90 | 2 +- src/physics/cam7/physpkg.F90 | 2 +- 4 files changed, 134 insertions(+), 153 deletions(-) diff --git a/src/physics/cam/beljaars_drag.F90 b/src/physics/cam/beljaars_drag.F90 index ccafd0e639..9d7fcfc6ef 100644 --- a/src/physics/cam/beljaars_drag.F90 +++ b/src/physics/cam/beljaars_drag.F90 @@ -1,152 +1,107 @@ +! Beljaars Sub-Grid Orographic (SGO) Form Drag Parameterization +! Returns drag profile and integrated stress associated with subgrid mountains +! with horizontal length scales nominally below 3 km. +! Based on Beljaars et al. (2004, QJRMS) https://doi.org/10.1256/qj.03.73. +! +! Original author: J. Bacmeister, March 2016, based on TMS module beljaars_drag - implicit none - private - save - - public init_blj ! Initialization - public compute_blj ! Full routine - - ! ------------ ! - ! Private data ! - ! ------------ ! - - integer, parameter :: r8 = selected_real_kind(12) ! 8 byte real + private - real(r8), parameter :: horomin= 1._r8 ! Minimum value of subgrid orographic height for mountain stress [ m ] - real(r8), parameter :: z0max = 100._r8 ! Maximum value of z_0 for orography [ m ] - real(r8), parameter :: dv2min = 0.01_r8 ! Minimum shear squared [ m2/s2 ] - real(r8) :: orocnst ! Converts from standard deviation to height [ no unit ] - real(r8) :: z0fac ! Factor determining z_0 from orographic standard deviation [ no unit ] - real(r8) :: karman ! von Karman constant - real(r8) :: gravit ! Acceleration due to gravity - real(r8) :: rair ! Gas constant for dry air + public :: beljaars_drag_run contains +!> \section arg_table_beljaars_drag_run Argument Table +!! \htmlinclude beljaars_drag_run.html + ! Compute Beljaars SGO form drag profile and surface stresses + subroutine beljaars_drag_run( & + do_beljaars, & + ncol, pver, & + u, v, t, pmid, delp, zm, sgh30, & + gravit, & + ! below output: + drag, taux, tauy, & + errmsg, errflg) + use ccpp_kinds, only: kind_phys + + logical, intent(in) :: do_beljaars ! is Beljaars active? + integer, intent(in) :: ncol + integer, intent(in) :: pver + real(kind_phys), intent(in) :: u(:, :) ! zonal wind [m s-1] + real(kind_phys), intent(in) :: v(:, :) ! meridional wind [m s-1] + real(kind_phys), intent(in) :: t(:, :) ! air temperature [K] + real(kind_phys), intent(in) :: pmid(:, :) ! air pressure [Pa] + real(kind_phys), intent(in) :: delp(:, :) ! air pressure thickness [Pa] + real(kind_phys), intent(in) :: zm(:, :) ! geopotential height wrt surface [m] + real(kind_phys), intent(in) :: sgh30(:) ! standard deviation of subgrid orography [m] + real(kind_phys), intent(in) :: gravit ! gravitational acceleration [m s-2] + + real(kind_phys), intent(out) :: drag(:, :) ! SGO drag profile [s-1] + real(kind_phys), intent(out) :: taux(:) ! surface zonal wind stress [N m-2] + real(kind_phys), intent(out) :: tauy(:) ! surface meridional wind stress [N m-2] + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Local variables + integer :: i, k + + real(kind_phys) :: vmag ! velocity magnitude [m s-1] + + real(kind_phys) :: alpha, beta, Cmd, Ccorr, n1, n2, k1, kflt, k2, IH + real(kind_phys) :: a1(ncol), a2(ncol) + + errmsg = '' + errflg = 0 + + if(.not. do_beljaars) then + ! if not doing Beljaars, return zero drag and stresses from routine. + drag(:,:) = 0._kind_phys + taux(:,:) = 0._kind_phys + tauy(:,:) = 0._kind_phys + return + end if + + alpha = 12.0_kind_phys + beta = 1.0_kind_phys + n1 = -1.9_kind_phys + n2 = -2.8_kind_phys + + Cmd = 0.005_kind_phys + Ccorr = 0.6_kind_phys * 5.0_kind_phys + + kflt = 0.00035_kind_phys ! m-1 + k1 = 0.003_kind_phys ! m-1 + IH = 0.00102_kind_phys ! m-1 + + a1(1:ncol) = (sgh30(1:ncol) * sgh30(1:ncol)) / (IH * (kflt**n1)) + a2(1:ncol) = a1(1:ncol) * k1**(n1 - n2) - !============================================================================ ! - ! ! - !============================================================================ ! - - subroutine init_blj( kind, gravit_in, rair_in , errstring ) - - integer, intent(in) :: kind - real(r8), intent(in) :: gravit_in, rair_in - - character(len=*), intent(out) :: errstring - - errstring = ' ' - - if ( kind /= r8 ) then - errstring = 'inconsistent KIND of reals passed to init_blj' - return - endif - - gravit = gravit_in - rair = rair_in - - end subroutine init_blj - - !============================================================================ ! - ! ! - !============================================================================ ! - - subroutine compute_blj( pcols , pver , ncol , & - u , v , t , pmid , delp , & - zm , sgh , drag , taux , tauy , & - landfrac ) - - !------------------------------------------------------------------------------ ! - ! Beljaars Sub-Grid Orographic (SGO) Form drag parameterization ! - ! ! - ! Returns drag profile and integrated stress associated with subgrid mountains ! - ! with horizontal length scales nominally below 3km. Similar to TMS but ! - ! drag is distributed in the vertical (Beljaars et al., 2003, QJRMS). ! - ! ! - ! First cut follows TMS. J. Bacmeister, March 2016 ! - !------------------------------------------------------------------------------ ! - - ! ---------------------- ! - ! Input-Output Arguments ! - ! ---------------------- ! - - integer, intent(in) :: pcols ! Number of columns dimensioned - integer, intent(in) :: pver ! Number of model layers - integer, intent(in) :: ncol ! Number of columns actually used - - real(r8), intent(in) :: u(pcols,pver) ! Layer mid-point zonal wind [ m/s ] - real(r8), intent(in) :: v(pcols,pver) ! Layer mid-point meridional wind [ m/s ] - real(r8), intent(in) :: t(pcols,pver) ! Layer mid-point temperature [ K ] - real(r8), intent(in) :: pmid(pcols,pver) ! Layer mid-point pressure [ Pa ] - real(r8), intent(in) :: delp(pcols,pver) ! Layer thickness [ Pa ] - real(r8), intent(in) :: zm(pcols,pver) ! Layer mid-point height [ m ] - real(r8), intent(in) :: sgh(pcols) ! Standard deviation of orography [ m ] - real(r8), intent(in) :: landfrac(pcols) ! Land fraction [ fraction ] - - real(r8), intent(out) :: drag(pcols,pver) ! SGO drag profile [ kg/s/m2 ] - real(r8), intent(out) :: taux(pcols) ! Surface zonal wind stress [ N/m2 ] - real(r8), intent(out) :: tauy(pcols) ! Surface meridional wind stress [ N/m2 ] - - ! --------------- ! - ! Local Variables ! - ! --------------- ! - - integer :: i,k ! Loop indices - integer :: kb, kt ! Bottom and top of source region - - real(r8) :: vmag ! Velocity magnitude [ m /s ] - - real(r8) :: alpha,beta,Cmd,Ccorr,n1,n2,k1,kflt,k2,IH - real(r8) :: a1(pcols),a2(pcols) - - alpha = 12._r8 - beta = 1._r8 - n1 = -1.9_r8 - n2 = -2.8_r8 - - Cmd = 0.005_r8 - Ccorr = 0.6_r8 * 5._r8 - - kflt = 0.00035_r8 ! m-1 - k1 = 0.003_r8 ! m-1 - IH = 0.00102_r8 ! m-1 - - a1(1:ncol) = (sgh(1:ncol)*sgh(1:ncol)) / ( IH* (kflt**n1) ) - a2(1:ncol) = a1(1:ncol) * k1**(n1-n2) - - - ! ----------------------- ! - ! Main Computation Begins ! - ! ----------------------- ! - do k = 1, pver - do i = 1, ncol - Vmag = SQRT( u(i,k)**2 + v(i,k)**2) - drag(i,k) = -alpha * beta * Cmd * Ccorr * Vmag * 2.109_r8 * & - EXP ( -(zm(i,k)/1500._r8 )*SQRT(zm(i,k)/1500._r8) ) * ( zm(i,k)**(-1.2_r8) ) & - * a2(i) + do i = 1, ncol + vmag = sqrt(u(i, k)**2 + v(i, k)**2) + drag(i, k) = -alpha * beta * Cmd * Ccorr * vmag * 2.109_kind_phys * & + exp(-(zm(i, k) / 1500.0_kind_phys) * sqrt(zm(i, k) / 1500.0_kind_phys)) * & + (zm(i, k)**(-1.2_kind_phys)) * a2(i) + end do end do - end do - - - !---------------------------------! - ! Diagnose effective surface drag ! - ! in X and Y by integrating in ! - ! the vertical ! - !---------------------------------! - ! FIXME: uses 'state' u and v. - ! Should updated u and v's be used? - - taux=0._r8 - tauy=0._r8 + + ! Diagnose effective surface drag in X and Y by integrating in the vertical + ! + ! taux, tauy stresses generated by Beljaars drag here + ! uses the pre-vertical diffusion winds (and not the updated winds) + ! however, at this point only these winds are available because the diffusion + ! solver runs after the orographic drag. + ! the actual provisionally updated winds are actually used to recompute taubljx/y + ! which are added to the updated residual stress. see the physics scheme + ! beljaars_add_updated_residual_stress. (hplin, 5/13/26) + taux(:) = 0.0_kind_phys + tauy(:) = 0.0_kind_phys do k = 1, pver - do i = 1, ncol - taux(i) = taux(i) + drag(i,k)*u(i,k)*delp(i,k)/gravit - tauy(i) = tauy(i) + drag(i,k)*v(i,k)*delp(i,k)/gravit + do i = 1, ncol + taux(i) = taux(i) + drag(i, k) * u(i, k) * delp(i, k) / gravit + tauy(i) = tauy(i) + drag(i, k) * v(i, k) * delp(i, k) / gravit + end do end do - end do - - return - end subroutine compute_blj + end subroutine beljaars_drag_run end module beljaars_drag diff --git a/src/physics/cam/beljaars_drag_cam.F90 b/src/physics/cam/beljaars_drag_cam.F90 index d81d2bb9b0..ab45050e4d 100644 --- a/src/physics/cam/beljaars_drag_cam.F90 +++ b/src/physics/cam/beljaars_drag_cam.F90 @@ -81,9 +81,7 @@ subroutine beljaars_drag_init() use cam_history, only: addfld, add_default, horiz_only use error_messages, only: handle_errmsg use phys_control, only: phys_getopts - use physconst, only: karman, gravit, rair use physics_buffer, only: pbuf_get_index - use beljaars_drag, only: init_blj logical :: history_amwg @@ -93,7 +91,6 @@ subroutine beljaars_drag_init() call phys_getopts(history_amwg_out=history_amwg) - call init_blj( r8, gravit, rair, errstring ) call handle_errmsg(errstring, subname="init_blj") call addfld('DRAGBLJ', (/ 'lev' /) , 'A', '1/s', 'Drag profile from Beljaars SGO ') @@ -112,21 +109,27 @@ subroutine beljaars_drag_init() end subroutine beljaars_drag_init -subroutine beljaars_drag_tend(state, pbuf, cam_in) +subroutine beljaars_drag_tend(state, pbuf) use physics_buffer, only: physics_buffer_desc, pbuf_get_field use physics_types, only: physics_state - use camsrfexch, only: cam_in_t use cam_history, only: outfld - use beljaars_drag, only: compute_blj + + use physconst, only: gravit + use beljaars_drag, only: beljaars_drag_run type(physics_state), intent(in) :: state type(physics_buffer_desc), pointer, intent(in) :: pbuf(:) - type(cam_in_t), intent(in) :: cam_in real(r8), pointer :: sgh30(:) real(r8), pointer :: dragblj(:,:) real(r8), pointer :: taubljx(:), taubljy(:) + integer :: ncol + character(len=512) :: errmsg + integer :: errflg + + ncol = state%ncol + call pbuf_get_field(pbuf, dragblj_idx, dragblj) call pbuf_get_field(pbuf, taubljx_idx, taubljx) call pbuf_get_field(pbuf, taubljy_idx, taubljy) @@ -140,10 +143,33 @@ subroutine beljaars_drag_tend(state, pbuf, cam_in) call pbuf_get_field(pbuf, sgh30_idx, sgh30) - call compute_blj( pcols , pver , state%ncol , & - state%u , state%v , state%t , state%pmid , & - state%pdel , state%zm , sgh30 , dragblj , & - taubljx , taubljy , cam_in%landfrac ) + ! zero to pcols + dragblj(:, :) = 0._r8 + taubljx(:) = 0._r8 + taubljy(:) = 0._r8 + + ! Call the CCPPized subroutine: + call beljaars_drag_run( & + do_beljaars = do_beljaars, & + ncol = state%ncol, & + pver = pver, & + u = state%u(:ncol, :), & + v = state%v(:ncol, :), & + t = state%t(:ncol, :), & + pmid = state%pmid(:ncol, :), & + delp = state%pdel(:ncol, :), & + zm = state%zm(:ncol, :), & + sgh30 = sgh30(:ncol), & + gravit = gravit, & + dragblj = dragblj(:ncol, :), & + taubljx = taubljx(:ncol), & + taubljy = taubljy(:ncol), & + errmsg = errmsg, & + errflg = errflg) + + if(errflg /= 0) then + call endrun('beljaars_drag_run: '//errmsg) + end if call outfld("TAUBLJX", taubljx, pcols, state%lchnk) call outfld("TAUBLJY", taubljy, pcols, state%lchnk) diff --git a/src/physics/cam/physpkg.F90 b/src/physics/cam/physpkg.F90 index 27790526e5..b4cc220f2e 100644 --- a/src/physics/cam/physpkg.F90 +++ b/src/physics/cam/physpkg.F90 @@ -1689,7 +1689,7 @@ subroutine tphysac (ztodt, cam_in, & end if call trb_mtn_stress_tend(state, pbuf, cam_in) - call beljaars_drag_tend(state, pbuf, cam_in) + call beljaars_drag_tend(state, pbuf) if (trim(cam_take_snapshot_after) == "orographic_form_drag_stress") then call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf,& diff --git a/src/physics/cam7/physpkg.F90 b/src/physics/cam7/physpkg.F90 index ae1620787a..a623c4c085 100644 --- a/src/physics/cam7/physpkg.F90 +++ b/src/physics/cam7/physpkg.F90 @@ -2204,7 +2204,7 @@ subroutine tphysac (ztodt, cam_in, & fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) end if - call beljaars_drag_tend(state, pbuf, cam_in) + call beljaars_drag_tend(state, pbuf) ! TMS is not active in CAM7 (it is only for CAM5), but the tms tend subroutine ! will initialize the pbuf fields to zero - no logic is computed below: From 0d43c6686a3eb2464c510fdda5d4ead79ba52f7f Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Wed, 13 May 2026 15:01:11 -0400 Subject: [PATCH 2/3] Move beljaars_drag.F90 to atmos_phys. --- bld/configure | 1 + src/physics/cam/beljaars_drag.F90 | 107 ------------------------------ 2 files changed, 1 insertion(+), 107 deletions(-) delete mode 100644 src/physics/cam/beljaars_drag.F90 diff --git a/bld/configure b/bld/configure index 6193e62b02..8a7a96dcc0 100755 --- a/bld/configure +++ b/bld/configure @@ -2150,6 +2150,7 @@ sub write_filepath print $fh "$camsrcdir/src/atmos_phys/schemes/cloud_fraction\n"; print $fh "$camsrcdir/src/atmos_phys/schemes/vertical_diffusion\n"; print $fh "$camsrcdir/src/atmos_phys/schemes/holtslag_boville\n"; + print $fh "$camsrcdir/src/atmos_phys/schemes/beljaars_drag\n"; # Dynamics package and test utilities print $fh "$camsrcdir/src/dynamics/$dyn\n"; diff --git a/src/physics/cam/beljaars_drag.F90 b/src/physics/cam/beljaars_drag.F90 deleted file mode 100644 index 9d7fcfc6ef..0000000000 --- a/src/physics/cam/beljaars_drag.F90 +++ /dev/null @@ -1,107 +0,0 @@ -! Beljaars Sub-Grid Orographic (SGO) Form Drag Parameterization -! Returns drag profile and integrated stress associated with subgrid mountains -! with horizontal length scales nominally below 3 km. -! Based on Beljaars et al. (2004, QJRMS) https://doi.org/10.1256/qj.03.73. -! -! Original author: J. Bacmeister, March 2016, based on TMS -module beljaars_drag - implicit none - private - - public :: beljaars_drag_run - -contains -!> \section arg_table_beljaars_drag_run Argument Table -!! \htmlinclude beljaars_drag_run.html - ! Compute Beljaars SGO form drag profile and surface stresses - subroutine beljaars_drag_run( & - do_beljaars, & - ncol, pver, & - u, v, t, pmid, delp, zm, sgh30, & - gravit, & - ! below output: - drag, taux, tauy, & - errmsg, errflg) - use ccpp_kinds, only: kind_phys - - logical, intent(in) :: do_beljaars ! is Beljaars active? - integer, intent(in) :: ncol - integer, intent(in) :: pver - real(kind_phys), intent(in) :: u(:, :) ! zonal wind [m s-1] - real(kind_phys), intent(in) :: v(:, :) ! meridional wind [m s-1] - real(kind_phys), intent(in) :: t(:, :) ! air temperature [K] - real(kind_phys), intent(in) :: pmid(:, :) ! air pressure [Pa] - real(kind_phys), intent(in) :: delp(:, :) ! air pressure thickness [Pa] - real(kind_phys), intent(in) :: zm(:, :) ! geopotential height wrt surface [m] - real(kind_phys), intent(in) :: sgh30(:) ! standard deviation of subgrid orography [m] - real(kind_phys), intent(in) :: gravit ! gravitational acceleration [m s-2] - - real(kind_phys), intent(out) :: drag(:, :) ! SGO drag profile [s-1] - real(kind_phys), intent(out) :: taux(:) ! surface zonal wind stress [N m-2] - real(kind_phys), intent(out) :: tauy(:) ! surface meridional wind stress [N m-2] - character(len=*), intent(out) :: errmsg - integer, intent(out) :: errflg - - ! Local variables - integer :: i, k - - real(kind_phys) :: vmag ! velocity magnitude [m s-1] - - real(kind_phys) :: alpha, beta, Cmd, Ccorr, n1, n2, k1, kflt, k2, IH - real(kind_phys) :: a1(ncol), a2(ncol) - - errmsg = '' - errflg = 0 - - if(.not. do_beljaars) then - ! if not doing Beljaars, return zero drag and stresses from routine. - drag(:,:) = 0._kind_phys - taux(:,:) = 0._kind_phys - tauy(:,:) = 0._kind_phys - return - end if - - alpha = 12.0_kind_phys - beta = 1.0_kind_phys - n1 = -1.9_kind_phys - n2 = -2.8_kind_phys - - Cmd = 0.005_kind_phys - Ccorr = 0.6_kind_phys * 5.0_kind_phys - - kflt = 0.00035_kind_phys ! m-1 - k1 = 0.003_kind_phys ! m-1 - IH = 0.00102_kind_phys ! m-1 - - a1(1:ncol) = (sgh30(1:ncol) * sgh30(1:ncol)) / (IH * (kflt**n1)) - a2(1:ncol) = a1(1:ncol) * k1**(n1 - n2) - - do k = 1, pver - do i = 1, ncol - vmag = sqrt(u(i, k)**2 + v(i, k)**2) - drag(i, k) = -alpha * beta * Cmd * Ccorr * vmag * 2.109_kind_phys * & - exp(-(zm(i, k) / 1500.0_kind_phys) * sqrt(zm(i, k) / 1500.0_kind_phys)) * & - (zm(i, k)**(-1.2_kind_phys)) * a2(i) - end do - end do - - ! Diagnose effective surface drag in X and Y by integrating in the vertical - ! - ! taux, tauy stresses generated by Beljaars drag here - ! uses the pre-vertical diffusion winds (and not the updated winds) - ! however, at this point only these winds are available because the diffusion - ! solver runs after the orographic drag. - ! the actual provisionally updated winds are actually used to recompute taubljx/y - ! which are added to the updated residual stress. see the physics scheme - ! beljaars_add_updated_residual_stress. (hplin, 5/13/26) - taux(:) = 0.0_kind_phys - tauy(:) = 0.0_kind_phys - do k = 1, pver - do i = 1, ncol - taux(i) = taux(i) + drag(i, k) * u(i, k) * delp(i, k) / gravit - tauy(i) = tauy(i) + drag(i, k) * v(i, k) * delp(i, k) / gravit - end do - end do - - end subroutine beljaars_drag_run -end module beljaars_drag From 2044264a0091b9c650f39ce179aaee2e3e7ce59d Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Wed, 13 May 2026 15:21:03 -0400 Subject: [PATCH 3/3] Remove unused arguments --- src/physics/cam/beljaars_drag_cam.F90 | 19 ++++--------------- 1 file changed, 4 insertions(+), 15 deletions(-) diff --git a/src/physics/cam/beljaars_drag_cam.F90 b/src/physics/cam/beljaars_drag_cam.F90 index ab45050e4d..6010d9d7ee 100644 --- a/src/physics/cam/beljaars_drag_cam.F90 +++ b/src/physics/cam/beljaars_drag_cam.F90 @@ -18,10 +18,6 @@ module beljaars_drag_cam ! Is this module on at all? logical, public, protected :: do_beljaars = .false. -! Tuning parameters for TMS. -real(r8) :: blj_orocnst -real(r8) :: blj_z0fac - ! pbuf field indices integer :: & sgh30_idx = -1, & @@ -34,7 +30,7 @@ module beljaars_drag_cam subroutine beljaars_drag_readnl(nlfile) use namelist_utils, only: find_group_name use units, only: getunit, freeunit - use spmd_utils, only: masterprocid, mpi_logical, mpi_real8, mpicom + use spmd_utils, only: masterprocid, mpi_logical, mpicom ! filepath for file containing namelist input character(len=*), intent(in) :: nlfile @@ -79,20 +75,15 @@ end subroutine beljaars_drag_register subroutine beljaars_drag_init() use cam_history, only: addfld, add_default, horiz_only - use error_messages, only: handle_errmsg use phys_control, only: phys_getopts use physics_buffer, only: pbuf_get_index logical :: history_amwg - character(len=128) :: errstring - if (.not. do_beljaars) return call phys_getopts(history_amwg_out=history_amwg) - call handle_errmsg(errstring, subname="init_blj") - call addfld('DRAGBLJ', (/ 'lev' /) , 'A', '1/s', 'Drag profile from Beljaars SGO ') call addfld('TAUBLJX', horiz_only, 'A', 'N/m2', 'Zonal integrated drag from Beljaars SGO') call addfld('TAUBLJY', horiz_only, 'A', 'N/m2', 'Meridional integrated drag from Beljaars SGO') @@ -155,15 +146,13 @@ subroutine beljaars_drag_tend(state, pbuf) pver = pver, & u = state%u(:ncol, :), & v = state%v(:ncol, :), & - t = state%t(:ncol, :), & - pmid = state%pmid(:ncol, :), & delp = state%pdel(:ncol, :), & zm = state%zm(:ncol, :), & sgh30 = sgh30(:ncol), & gravit = gravit, & - dragblj = dragblj(:ncol, :), & - taubljx = taubljx(:ncol), & - taubljy = taubljy(:ncol), & + drag = dragblj(:ncol, :), & + taux = taubljx(:ncol), & + tauy = taubljy(:ncol), & errmsg = errmsg, & errflg = errflg)