Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

implementation of new orographic drag parameterization scheme to alleviate wind bias in E3SM #6665

Closed
wants to merge 2 commits into from
Closed

implementation of new orographic drag parameterization scheme to alleviate wind bias in E3SM #6665

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
74a116f
Subject (<= 72 Characters)
jinboxie Oct 2, 2024
392d1d9
Subject (<= 72 Characters)
jinboxie Oct 2, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion components/eam/bld/namelist_files/namelist_defaults_eam.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -133,7 +133,7 @@
<bnd_topo hgrid="ne16np4" npg="0">atm/cam/topo/USGS-gtopo30_ne16np4_16xconsistentSGH.c20160612.nc</bnd_topo>
<bnd_topo hgrid="ne16np4" npg="2">atm/cam/topo/USGS-gtopo30_ne16np4pg2_16xdel2_20200527.nc</bnd_topo>
<bnd_topo hgrid="ne30np4" npg="0">atm/cam/topo/USGS-gtopo30_ne30np4_16xdel2-PFC-consistentSGH.nc</bnd_topo>
<bnd_topo hgrid="ne30np4" npg="2">atm/cam/topo/USGS-gtopo30_ne30np4pg2_x6t-SGH.c20210614.nc</bnd_topo>
<bnd_topo hgrid="ne30np4" npg="2">atm/cam/topo/USGS-gtopo30_ne30np4pg2_x6t-SGH_forOroDrag.c20241001.nc</bnd_topo>
<bnd_topo hgrid="ne30np4" npg="3">atm/cam/topo/USGS-gtopo30_ne30np4pg3_16xdel2.c20200504.nc</bnd_topo>
<bnd_topo hgrid="ne30np4" npg="4">atm/cam/topo/USGS-gtopo30_ne30np4pg4_16xdel2.c20200504.nc</bnd_topo>
<bnd_topo hgrid="ne45np4" npg="2">atm/cam/topo/USGS-gtopo30_ne45np4pg2_16xdel2.c20200615.nc</bnd_topo>
Expand Down
41 changes: 41 additions & 0 deletions components/eam/src/control/startup_initialconds.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,16 +5,28 @@ module startup_initialconds
!
!-----------------------------------------------------------------------

use pio, only: file_desc_t

implicit none
private
save

public :: initial_conds ! Read in initial conditions (dycore dependent)
!added for orographic drag
public topoGWD_file_get_id
public setup_initialGWD
public close_initial_fileGWD
type(file_desc_t), pointer :: ncid_topoGWD

!=======================================================================
contains
!=======================================================================

function topoGWD_file_get_id()
type(file_desc_t), pointer :: topoGWD_file_get_id
topoGWD_file_get_id => ncid_topoGWD
end function topoGWD_file_get_id

subroutine initial_conds(dyn_in)

! This routine does some initializing of buffers that should move to a
Expand Down Expand Up @@ -62,4 +74,33 @@ end subroutine initial_conds

!=======================================================================

subroutine setup_initialGWD()
use filenames, only: bnd_topo
use ioFileMod, only: getfil
use cam_pio_utils, only: cam_pio_openfile
use pio, only: pio_nowrite
!
! Input arguments
!
!-----------------------------------------------------------------------
include 'netcdf.inc'
!-----------------------------------------------------------------------
character(len=256) :: bnd_topo_loc ! filepath of topo file on local disk
allocate(ncid_topoGWD)
call getfil(bnd_topo, bnd_topo_loc)
call cam_pio_openfile(ncid_topoGWD, bnd_topo_loc, PIO_NOWRITE)
end subroutine setup_initialGWD

subroutine close_initial_fileGWD
use pio, only: pio_closefile
call pio_closefile(ncid_topoGWD)
deallocate(ncid_topoGWD)
nullify(ncid_topoGWD)
end subroutine close_initial_fileGWD
!=======================================================================





end module startup_initialconds
125 changes: 119 additions & 6 deletions components/eam/src/physics/cam/clubb_intr.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -927,7 +927,18 @@ subroutine clubb_ini_cam(pbuf2d, dp1_in)
call addfld ('VMAGDP', horiz_only, 'A', '-', 'ZM gustiness enhancement')
call addfld ('VMAGCL', horiz_only, 'A', '-', 'CLUBB gustiness enhancement')
call addfld ('TPERTBLT', horiz_only, 'A', 'K', 'perturbation temperature at PBL top')

!==================================
!!added for TOFD output
call addfld ('DTAUX3_FD',(/'lev'/),'A','m/s2','U tendency - fd orographic drag')
call addfld ('DTAUY3_FD',(/'lev'/),'A','m/s2','V tendency - fd orographic drag')
call addfld ('DUSFC_FD',horiz_only,'A','N/m2','fd zonal oro surface stress')
call addfld ('DVSFC_FD',horiz_only,'A','N/m2','fd merio oro surface stress')
call add_default('DTAUX3_FD', 1, ' ')
call add_default('DTAUY3_FD', 1, ' ')
call add_default('DUSFC_FD', 1, ' ')
call add_default('DVSFC_FD', 1, ' ')
!!added for TOFD output
!=====================================
! Initialize statistics, below are dummy variables
dum1 = 300._r8
dum2 = 1200._r8
Expand Down Expand Up @@ -1155,7 +1166,11 @@ subroutine clubb_tend_cam( &
use model_flags, only: ipdf_call_placement
use advance_clubb_core_module, only: ipdf_post_advance_fields
#endif

use gw_common, only: gwdo_gsd,grid_size,pblh_get_level_idx
use hycoef, only: etamid
use physconst, only: rh2o,pi,rearth,r_universal
!!get the znu,znw,p_top set to 0
use phys_grid, only: get_rlat_all_p
implicit none

! --------------- !
Expand Down Expand Up @@ -1518,7 +1533,24 @@ subroutine clubb_tend_cam( &

real(r8) :: sfc_v_diff_tau(pcols) ! Response to tau perturbation, m/s
real(r8), parameter :: pert_tau = 0.1_r8 ! tau perturbation, Pa

!===========================
!simply add par
!for z,dz,from other files
real(r8) :: ztop(pcols,pver) ! top interface height asl(m)
real(r8) :: zbot(pcols,pver) ! bottom interface height asl(m)
real(r8) :: zmid(pcols,pver) ! middle interface height asl(m)
real(r8) :: dz(pcols,pver)
real(r8) :: rlat(pcols) ! latitude in radians for columns
integer :: kpbl2d_in(pcols)
real(r8) :: ttgw(pcols,pver) ! temperature tendency
real(r8) :: utgw(pcols,pver) ! zonal wind tendency
real(r8) :: vtgw(pcols,pver) ! meridional wind tendency
real(r8) :: dtaux3_fd(pcols,pver)
real(r8) :: dtauy3_fd(pcols,pver)
real(r8) :: dusfc_fd(pcols)
real(r8) :: dvsfc_fd(pcols)
real(r8) :: dx(pcols),dy(pcols)
!==============================

real(r8) :: inv_exner_clubb_surf

Expand Down Expand Up @@ -1946,7 +1978,73 @@ subroutine clubb_tend_cam( &
tautmsx, tautmsy, cam_in%landfrac )
call t_stopf('compute_tms')
endif

ztop= 0.0_r8 ! top interface height asl(m)
zbot= 0.0_r8 ! bottom interface height asl(m)
zmid= 0.0_r8 ! middle interface height asl(m)
dz= 0.0_r8
kpbl2d_in = -1
dtaux3_fd= 0.0_r8
dtauy3_fd= 0.0_r8
dusfc_fd= 0.0_r8
dvsfc_fd= 0.0_r8
!similar as in gw_drag
do k=1,pverp-1
! assign values from top
ztop(1:ncol,k)=state%zi(1:ncol,pverp-k)
! assign values from bottom
zbot(1:ncol,k)=state%zi(1:ncol,pverp-k+1)
end do
!transform adding the pressure
!transfer from surface to sea level
do k=1,pver
do i=1,ncol
ztop(i,k)=ztop(i,k)+state%phis(i)/gravit
zbot(i,k)=zbot(i,k)+state%phis(i)/gravit
zmid(i,k)=state%zm(i,k)+state%phis(i)/gravit
!dz is from bottom to top already for gw_drag
dz(i,k)=ztop(i,k)-zbot(i,k)
end do
end do
!get the layer index of pblh in layer
kpbl2d_in=0._r8
do i=1,pcols
kpbl2d_in(i)=pblh_get_level_idx(zbot(i,:)-(state%phis(i)/gravit),pblh(i))
end do
!rlat
call get_rlat_all_p(lchnk, ncol, rlat)
!=========================================
utgw=0._r8
vtgw=0._r8
ttgw=0._r8
dusfc_fd=0._r8
dvsfc_fd=0._r8
!
call grid_size(state,dx,dy)
call gwdo_gsd(&
u3d=state%u(:,pver:1:-1),v3d=state%v(:,pver:1:-1),&
t3d=state%t(:,pver:1:-1),qv3d=state%q(:,pver:1:-1,1),&
p3d=state%pmid(:,pver:1:-1),p3di=state%pint(:,pver+1:1:-1),&
pi3d=state%exner(:,pver:1:-1),z=zbot,&
rublten=utgw(:,pver:1:-1),rvblten=vtgw(:,pver:1:-1),&
rthblten=ttgw(:,pver:1:-1),&
dtaux3d_fd=dtaux3_fd(:,pver:1:-1),dtauy3d_fd=dtauy3_fd(:,pver:1:-1),&
dusfcg_fd=dusfc_fd(:ncol),dvsfcg_fd=dvsfc_fd(:ncol),&
xland=cam_in%landfrac,br=state%ribulk,&
var2d=sgh30(:ncol),&
znu=etamid(pver:1:-1),dz=dz,pblh=pblh,&
cp=cpair,g=gravit,rd=rair,rv=rh2o,ep1=zvir,pi=pi,&
dx=dx,dy=dy,&
kpbl2d=kpbl2d_in,itimestep=hdtime,gwd_opt=0,&
ids=1,ide=pcols,jds=0,jde=0,kds=1,kde=pver, &
ims=1,ime=pcols,jms=0,jme=0,kms=1,kme=pver, &
its=1,ite=pcols,jts=0,jte=0,kts=1,kte=pver,&
gwd_ls=0,gwd_bl=0,gwd_ss=0,gwd_fd=1)
!!
call outfld ('DTAUX3_FD', dtaux3_fd, pcols, lchnk)
call outfld ('DTAUY3_FD', dtauy3_fd, pcols, lchnk)
call outfld ('DUSFC_FD', dusfc_fd, pcols, lchnk)
call outfld ('DVSFC_FD', dvsfc_fd, pcols, lchnk)
!!
if (micro_do_icesupersat) then
call physics_ptend_init(ptend_loc,state%psetcols, 'clubb_ice3', ls=.true., lu=.true., lv=.true., lq=lq)
endif
Expand Down Expand Up @@ -2067,7 +2165,12 @@ subroutine clubb_tend_cam( &
dum_core_rknd = real((ksrftms(i)*state1%v(i,pver)), kind = core_rknd)
vpwp_sfc = vpwp_sfc-(dum_core_rknd/rho_ds_zm(1))
endif

!----------------------------------------------------!
!Apply TOFD
!----------------------------------------------------!
!tendency is flipped already
um_forcing(2:pverp)=dtaux3_fd(i,pver:1:-1)
vm_forcing(2:pverp)=dtauy3_fd(i,pver:1:-1)
! Need to flip arrays around for CLUBB core
do k=1,pverp
um_in(k) = real(um(i,pverp-k+1), kind = core_rknd)
Expand Down Expand Up @@ -3112,6 +3215,7 @@ subroutine clubb_surface (state, cam_in, ustar, obklen)
use ppgrid, only: pver, pcols
use constituents, only: cnst_get_ind
use camsrfexch, only: cam_in_t
use hb_diff, only: pblintd_ri

implicit none

Expand Down Expand Up @@ -3143,6 +3247,7 @@ subroutine clubb_surface (state, cam_in, ustar, obklen)
real(r8) :: kinheat ! kinematic surface heat flux
real(r8) :: kinwat ! kinematic surface vapor flux
real(r8) :: kbfs ! kinematic surface buoyancy flux
real(r8) :: kbfs_pcol(pcols)
integer :: ixq,ixcldliq !PMA fix for thv
real(r8) :: rrho ! Inverse air density

Expand Down Expand Up @@ -3180,7 +3285,15 @@ subroutine clubb_surface (state, cam_in, ustar, obklen)
call calc_obklen( th(i), thv(i), cam_in%cflx(i,1), cam_in%shf(i), rrho, ustar(i), &
kinheat, kinwat, kbfs, obklen(i) )
enddo

!!===== add calculation of ribulk here=====
kbfs_pcol=0.0_r8
do i=1,ncol
call calc_obklen( th(i), thv(i), cam_in%cflx(i,1), cam_in%shf(i), rrho, ustar(i), &
kinheat, kinwat, kbfs, obklen(i) )
kbfs_pcol(i)=kbfs
enddo
call pblintd_ri(ncol, thv, state%zm, state%u, state%v, &
ustar, obklen, kbfs_pcol, state%ribulk)
return

#endif
Expand Down
41 changes: 39 additions & 2 deletions components/eam/src/physics/cam/comsrf.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ module comsrf
! USES:
!
use shr_kind_mod, only: r8 => shr_kind_r8, r4 => shr_kind_r4
use ppgrid, only: pcols, begchunk, endchunk
use ppgrid, only: pcols, begchunk, endchunk,nvar_dirOA,nvar_dirOL,indexb
use infnan, only: nan, assignment(=)
use cam_abortutils, only: endrun

Expand All @@ -31,6 +31,8 @@ module comsrf
! ! PUBLIC MEMBER FUNCTIONS:
!
public initialize_comsrf ! Set the surface temperature and sea-ice fraction
!!added for separate input of ogwd parareters in gw_drag
public initialize_comsrf2
!
! Public data
!
Expand All @@ -53,7 +55,14 @@ module comsrf
real(r8), allocatable:: prcsnw(:,:) ! cam tot snow precip
real(r8), allocatable:: trefmxav(:,:) ! diagnostic: tref max over the day
real(r8), allocatable:: trefmnav(:,:) ! diagnostic: tref min over the day

!!
public var,var30,oc,ol,oadir
real(r8), allocatable:: var(:,:)!sgh
real(r8), allocatable:: var30(:,:)!sgh30
real(r8), allocatable:: oc(:,:) ! Convexity
real(r8), allocatable:: oadir(:,:,:) ! Asymmetry
real(r8), allocatable:: ol(:,:,:) ! Effective length
!!
! Private module data

!===============================================================================
Expand Down Expand Up @@ -134,4 +143,32 @@ subroutine initialize_comsrf
end if
end subroutine initialize_comsrf

subroutine initialize_comsrf2
use cam_control_mod, only: ideal_phys, adiabatic
!-----------------------------------------------------------------------
!
! Purpose:
! Initialize surface data
!
! Method:
!
! Author: Mariana Vertenstein
!
!-----------------------------------------------------------------------
integer k,c ! level, constituent indices

if(.not. (adiabatic .or. ideal_phys)) then
allocate (var(pcols,begchunk:endchunk))
allocate (var30(pcols,begchunk:endchunk))
allocate (oc(pcols,begchunk:endchunk))
allocate (oadir(pcols,nvar_dirOA,begchunk:endchunk))
allocate (ol(pcols,nvar_dirOL,begchunk:endchunk))
var(:,:)=nan
var30(:,:)=nan
oc (:,:) = nan
oadir (:,:,:) = nan
ol (:,:,:) = nan
end if
end subroutine initialize_comsrf2

end module comsrf
Loading
Loading