diff --git a/components/eam/src/control/startup_initialconds.F90 b/components/eam/src/control/startup_initialconds.F90 index a68195c731db..68f9a2f12a3b 100644 --- a/components/eam/src/control/startup_initialconds.F90 +++ b/components/eam/src/control/startup_initialconds.F90 @@ -5,28 +5,16 @@ 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 topo_OD_file_get_id -public setup_initial_OD -public close_initial_file_OD -type(file_desc_t), pointer :: ncid_topo_OD !======================================================================= contains !======================================================================= -function topo_OD_file_get_id() - type(file_desc_t), pointer :: topo_OD_file_get_id - topo_OD_file_get_id => ncid_topo_OD -end function topo_OD_file_get_id - subroutine initial_conds(dyn_in) ! This routine does some initializing of buffers that should move to a @@ -72,35 +60,4 @@ subroutine initial_conds(dyn_in) end subroutine initial_conds -!======================================================================= - -subroutine setup_initial_OD() - 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_topo_OD) - call getfil(bnd_topo, bnd_topo_loc) - call cam_pio_openfile(ncid_topo_OD, bnd_topo_loc, PIO_NOWRITE) -end subroutine setup_initial_OD - -subroutine close_initial_file_OD - use pio, only: pio_closefile - call pio_closefile(ncid_topo_OD) - deallocate(ncid_topo_OD) - nullify(ncid_topo_OD) -end subroutine close_initial_file_OD -!======================================================================= - - - - - end module startup_initialconds diff --git a/components/eam/src/physics/cam/clubb_intr.F90 b/components/eam/src/physics/cam/clubb_intr.F90 index c69d76def44d..0d45232f8a98 100644 --- a/components/eam/src/physics/cam/clubb_intr.F90 +++ b/components/eam/src/physics/cam/clubb_intr.F90 @@ -2139,14 +2139,14 @@ 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 - if (use_od_fd) then + ! ------------------------------------------------- ! + ! Apply TOFD + ! ------------------------------------------------- ! + ! tendency is flipped already + if (use_od_fd) then um_forcing(2:pverp)=dtaux3_fd(i,pver:1:-1) vm_forcing(2:pverp)=dtauy3_fd(i,pver:1:-1) - endif + endif ! Need to flip arrays around for CLUBB core do k=1,pverp um_in(k) = real(um(i,pverp-k+1), kind = core_rknd) @@ -3170,7 +3170,7 @@ end subroutine clubb_tend_cam ! ! ! =============================================================================== ! - subroutine clubb_surface (state, cam_in, ustar, obklen) + subroutine clubb_surface (state, cam_in, pbuf, ustar, obklen) !------------------------------------------------------------------------------- ! Description: Provide the obukhov length and the surface friction velocity @@ -3192,7 +3192,7 @@ subroutine clubb_surface (state, cam_in, ustar, obklen) use constituents, only: cnst_get_ind use camsrfexch, only: cam_in_t use hb_diff, only: pblintd_ri - + use physics_buffer, only: pbuf_get_index, pbuf_get_field, physics_buffer_desc implicit none @@ -3200,8 +3200,9 @@ subroutine clubb_surface (state, cam_in, ustar, obklen) ! Input Auguments ! ! --------------- ! - type(physics_state), intent(inout) :: state ! Physics state variables - type(cam_in_t), intent(in) :: cam_in + type(physics_state), intent(inout) :: state ! Physics state variables + type(cam_in_t), intent(in) :: cam_in + type(physics_buffer_desc), pointer, intent(in) :: pbuf(:) ! ---------------- ! ! Output Auguments ! @@ -3231,6 +3232,9 @@ subroutine clubb_surface (state, cam_in, ustar, obklen) integer :: ixq,ixcldliq !PMA fix for thv real(r8) :: rrho ! Inverse air density + integer :: oro_drag_ribulk_idx ! pbuf index of bulk richardson number for oro drag + real(r8), pointer :: oro_drag_ribulk(:) ! pbuf pointer for bulk richardson number + #endif obklen(pcols) = 0.0_r8 @@ -3295,9 +3299,12 @@ subroutine clubb_surface (state, cam_in, ustar, obklen) kbfs_pcol(i)=kbfs enddo + oro_drag_ribulk_idx = pbuf_get_index('oro_drag_ribulk') + call pbuf_get_field(pbuf, oro_drag_ribulk_idx, oro_drag_ribulk) + !calculate the bulk richardson number call pblintd_ri(ncol, gravit, thv_lv, state%zm, state%u, state%v, & - ustar, obklen, kbfs_pcol, state%ribulk) + ustar, obklen, kbfs_pcol, oro_drag_ribulk) endif return diff --git a/components/eam/src/physics/cam/comsrf.F90 b/components/eam/src/physics/cam/comsrf.F90 index 64e3750dd4e7..7ac806c10326 100644 --- a/components/eam/src/physics/cam/comsrf.F90 +++ b/components/eam/src/physics/cam/comsrf.F90 @@ -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,nvar_dirOA,nvar_dirOL + use ppgrid, only: pcols, begchunk, endchunk use infnan, only: nan, assignment(=) use cam_abortutils, only: endrun @@ -31,8 +31,6 @@ 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_comsrf_OD ! ! Public data ! @@ -56,10 +54,6 @@ module comsrf real(r8), allocatable:: trefmxav(:,:) ! diagnostic: tref max over the day real(r8), allocatable:: trefmnav(:,:) ! diagnostic: tref min over the day - public oc, ol, oadir - real(r8), allocatable:: oc(:,:) ! Convexity - real(r8), allocatable:: oadir(:,:,:) ! Asymmetry - real(r8), allocatable:: ol(:,:,:) ! Effective length ! ! Private module data @@ -138,28 +132,4 @@ subroutine initialize_comsrf end if end subroutine initialize_comsrf - subroutine initialize_comsrf_OD - 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 (oc (pcols,begchunk:endchunk)) - allocate (oadir (pcols,nvar_dirOA,begchunk:endchunk)) - allocate (ol (pcols,nvar_dirOL,begchunk:endchunk)) - oc (:,:) = nan - oadir (:,:,:) = nan - ol (:,:,:) = nan - end if - end subroutine initialize_comsrf_OD - end module comsrf diff --git a/components/eam/src/physics/cam/gw_drag.F90 b/components/eam/src/physics/cam/gw_drag.F90 index 352858905ba9..96ac9f70021e 100644 --- a/components/eam/src/physics/cam/gw_drag.F90 +++ b/components/eam/src/physics/cam/gw_drag.F90 @@ -24,7 +24,7 @@ module gw_drag !-------------------------------------------------------------------------- use shr_kind_mod, only: r8 => shr_kind_r8 - use ppgrid, only: pcols, pver, pverp, nvar_dirOA, nvar_dirOL, begchunk, endchunk + use ppgrid, only: pcols, pver use hycoef, only: hyai, hybi, hyam, hybm, etamid use constituents, only: pcnst use physics_types, only: physics_state, physics_ptend, physics_ptend_init @@ -49,6 +49,7 @@ module gw_drag ! PUBLIC: interfaces ! public :: gw_drag_readnl ! Read namelist + public :: gw_register ! Register pbuf variables public :: gw_init ! Initialization public :: gw_tend ! interface to actual parameterization @@ -199,7 +200,16 @@ end subroutine gw_drag_readnl !========================================================================== -subroutine gw_init() +subroutine gw_register() + use od_common, only: oro_drag_register + + call oro_drag_register() + +end subroutine gw_register + +!========================================================================== + +subroutine gw_init(pbuf2d) !----------------------------------------------------------------------- ! Time independent initialization for multiple gravity wave ! parameterization. @@ -208,7 +218,7 @@ subroutine gw_init() use cam_history, only: addfld, horiz_only, add_default use interpolate_data, only: lininterp use phys_control, only: phys_getopts - use physics_buffer, only: pbuf_get_index + use physics_buffer, only: pbuf_get_index, physics_buffer_desc use ref_pres, only: pref_edge use physconst, only: gravit, rair @@ -218,12 +228,9 @@ subroutine gw_init() use gw_front, only: gw_front_init use gw_convect, only: gw_convect_init - use comsrf, only: oc, oadir, ol, initialize_comsrf_OD - use pio, only: file_desc_t - use startup_initialconds,only: topo_OD_file_get_id, setup_initial_OD, close_initial_file_OD - use ncdio_atm, only: infld - use cam_grid_support, only: cam_grid_check, cam_grid_get_decomp, cam_grid_id,cam_grid_get_dim_names - + use od_common, only: oro_drag_init + !------------------------------Arguments-------------------------------- + type(physics_buffer_desc), pointer :: pbuf2d(:,:) !---------------------------Local storage------------------------------- integer :: l, k @@ -296,36 +303,8 @@ subroutine gw_init() character(len=128) :: errstring !----------------------------------------------------------------------- - !added for input of od parameters - type(file_desc_t), pointer :: ncid_topo_OD - logical :: found=.false. - character(len=8) :: dim1name, dim2name - character*11 :: subname='gw_init' ! subroutine name - integer :: grid_id - pblh_idx = pbuf_get_index('pblh') - grid_id = cam_grid_id('physgrid') - - if (use_od_ls.or.use_od_bl) then - if (.not. cam_grid_check(grid_id)) then - call endrun(trim(subname)//': Internal error, no "physgrid" grid') - end if - call cam_grid_get_dim_names(grid_id, dim1name, dim2name) - call initialize_comsrf_OD() - call setup_initial_OD() - - ncid_topo_OD=>topo_OD_file_get_id() - call infld('OC', ncid_topo_OD, dim1name, dim2name, 1, pcols, begchunk, & - endchunk, oc , found, gridname='physgrid') - !keep the same interval of OA,OL - call infld('OA', ncid_topo_OD,dim1name, 'nvar_dirOA', dim2name, 1, pcols, 1, nvar_dirOA, begchunk, & - endchunk, oadir(:,:,:), found, gridname='physgrid') - call infld('OL', ncid_topo_OD,dim1name, 'nvar_dirOL', dim2name, 1, pcols, 1, nvar_dirOL, begchunk, & - endchunk, ol , found, gridname='physgrid') - if(.not. found) call endrun('ERROR: OD topo file readerr') - call close_initial_file_OD() - - endif + call oro_drag_init(pbuf2d) ! Set model flags. do_spectral_waves = (pgwv > 0 .and. (use_gw_front .or. use_gw_convect)) @@ -699,9 +678,6 @@ subroutine gw_tend(state, sgh, pbuf, dt, ptend, cam_in) ! real(r8), pointer :: pblh(:) real(r8) :: dx(pcols),dy(pcols) - ! - logical :: gwd_ls,gwd_bl,gwd_ss,gwd_fd - ! !---------------------------Local storage------------------------------- @@ -998,22 +974,12 @@ subroutine gw_tend(state, sgh, pbuf, dt, ptend, cam_in) ttgw, qtgw, taucd, egwdffi, gwut(:,:,0:0), dttdf, dttke) endif ! - if (use_od_ls.or.& - use_od_bl.or.& - use_od_ss) then - !open ogwd,bl,ss, - !close fd - gwd_ls=use_od_ls - gwd_bl=use_od_bl - gwd_ss=use_od_ss - gwd_fd=.false. - ! + if ( use_od_ls .or. use_od_bl .or. use_od_ss) then utgw=0.0_r8 vtgw=0.0_r8 ttgw=0.0_r8 - ! call oro_drag_interface(state,cam_in,sgh,pbuf,dt,nm,& - gwd_ls,gwd_bl,gwd_ss,gwd_fd,& + use_od_ls,use_od_bl,use_od_ss,.false.,& od_ls_ncleff,od_bl_ncd,od_ss_sncleff,& utgw,vtgw,ttgw,& dtaux3_ls=dtaux3_ls,dtauy3_ls=dtauy3_ls,& @@ -1024,14 +990,13 @@ subroutine gw_tend(state, sgh, pbuf, dt, ptend, cam_in) dusfc_bl=dusfc_bl,dvsfc_bl=dvsfc_bl,& dusfc_ss=dusfc_ss,dvsfc_ss=dvsfc_ss,& dusfc_fd=dummx_fd,dvsfc_fd=dummy_fd) - endif - ! - ! Add the orographic tendencies to the spectrum tendencies - ! Compute the temperature tendency from energy conservation - ! (includes spectrum). - ! both old and new gwd scheme will add the tendency to circulation - ! + ! + ! Add the orographic tendencies to the spectrum tendencies + ! Compute the temperature tendency from energy conservation + ! (includes spectrum). + ! both old and new gwd scheme will add the tendency to circulation + ! if (use_gw_oro.or.& use_od_ls .or.& use_od_bl .or.& diff --git a/components/eam/src/physics/cam/od_common.F90 b/components/eam/src/physics/cam/od_common.F90 index 3eb81889e95a..5d84e718b287 100644 --- a/components/eam/src/physics/cam/od_common.F90 +++ b/components/eam/src/physics/cam/od_common.F90 @@ -1,5 +1,5 @@ module od_common -! +!========================================================================== ! This module contains code common to different orographic drag ! parameterizations. ! It includes 4 parts: @@ -7,23 +7,159 @@ module od_common ! flow-blocking drag (Xie et al.,2020), ! small-scale orographic gravity wave drag (Tsiringakis et al. 2017), ! turbulent orographic form drag (Beljaars et al.,2004). -! -use gw_utils, only: r8 -use ppgrid, only: pver,nvar_dirOA,nvar_dirOL +!========================================================================== +use shr_kind_mod, only: i8 => shr_kind_i8, r8 => shr_kind_r8 +use shr_sys_mod, only: shr_sys_flush +use ppgrid, only: pcols, pver, begchunk, endchunk use cam_logfile, only: iulog +use cam_abortutils,only: endrun +use pio, only: file_desc_t +use phys_control, only: use_od_ls, use_od_bl, use_od_ss, od_ls_ncleff, od_bl_ncd, od_ss_sncleff +use physics_buffer,only: dtype_r8, physics_buffer_desc, pbuf_get_chunk +use physics_buffer,only: pbuf_get_index, pbuf_get_field, pbuf_add_field, pbuf_set_field implicit none private save ! Public interface. +public :: oro_drag_register +public :: oro_drag_init public :: oro_drag_interface public :: od_gsd,pblh_get_level_idx,grid_size +type(file_desc_t), pointer :: topo_file_ncid + +! dimensions for topo shape data +integer, parameter :: ndir_asymmetry = 2+1 ! add 1 to avoid bug reading file - not sure why this happens +integer, parameter :: ndir_efflength = 180 ! 1-degree resolution with opposite directions mirrored + +! pbuf indices for data read in from topo data file +integer :: oro_drag_convexity_idx = -1 ! Convexity +integer :: oro_drag_asymmetry_idx = -1 ! Asymmetry +integer :: oro_drag_efflength_idx = -1 ! Effective length +integer :: oro_drag_ribulk_idx = -1 ! bulk richardson number (calculated in CLUBB) + contains !========================================================================== +subroutine oro_drag_open_topo_file() + use filenames, only: bnd_topo + use ioFileMod, only: getfil + use cam_pio_utils,only: cam_pio_openfile + use pio, only: pio_nowrite + include 'netcdf.inc' + !----------------------------------------------------------------------- + character(len=256) :: bnd_topo_loc ! filepath of topo file on local disk + allocate(topo_file_ncid) + call getfil(bnd_topo, bnd_topo_loc) + call cam_pio_openfile(topo_file_ncid, bnd_topo_loc, PIO_NOWRITE) +end subroutine oro_drag_open_topo_file + +!========================================================================== + +subroutine oro_drag_close_topo_file + use pio, only: pio_closefile + call pio_closefile(topo_file_ncid) + deallocate(topo_file_ncid) + nullify(topo_file_ncid) +end subroutine oro_drag_close_topo_file + +!========================================================================== + +subroutine oro_drag_register() + !----------------------------------------------------------------------- + ! Register pbuf variables for orographic drag parameterizations + !----------------------------------------------------------------------- + ! create pbuf variables to hold oro drag data + if (use_od_ls.or.use_od_bl) then + call pbuf_add_field('oro_drag_convexity','physpkg',dtype_r8,(/pcols/), oro_drag_convexity_idx) + call pbuf_add_field('oro_drag_asymmetry','physpkg',dtype_r8,(/pcols,ndir_asymmetry/),oro_drag_asymmetry_idx) + call pbuf_add_field('oro_drag_efflength','physpkg',dtype_r8,(/pcols,ndir_efflength/),oro_drag_efflength_idx) + end if + if (use_od_ss) then + call pbuf_add_field('oro_drag_ribulk', 'physpkg',dtype_r8,(/pcols/), oro_drag_ribulk_idx) + end if + +end subroutine oro_drag_register + +!========================================================================== + +subroutine oro_drag_init(pbuf2d) + !----------------------------------------------------------------------- + ! Initialization for orographic drag parameterizations + !----------------------------------------------------------------------- + use pio, only: file_desc_t + use ncdio_atm, only: infld + use cam_grid_support, only: cam_grid_check, cam_grid_get_decomp, cam_grid_id,cam_grid_get_dim_names + use infnan, only: nan, assignment(=) + !----------------------------------------------------------------------- + type(physics_buffer_desc), pointer :: pbuf2d(:,:) + !----------------------------------------------------------------------- + logical :: found + character(len=8) :: dim1name, dim2name + character*11 :: subname='oro_drag_init' + integer :: grid_id + integer :: c + + real(r8), allocatable:: oro_drag_convexity_tmp(:,:) + real(r8), allocatable:: oro_drag_asymmetry_tmp(:,:,:) + real(r8), allocatable:: oro_drag_efflength_tmp(:,:,:) + + type(physics_buffer_desc), pointer :: pbuf_chunk(:) ! temporary pbuf pointer for single chunk + !----------------------------------------------------------------------- + if (.not.(use_od_ls.or.use_od_bl)) return + + grid_id = cam_grid_id('physgrid') + if (.not. cam_grid_check(grid_id)) then + call endrun(trim(subname)//': Internal error, no "physgrid" grid') + end if + + ! Alocate variables for reading oro drag data + allocate( oro_drag_convexity_tmp(pcols,begchunk:endchunk) ) + allocate( oro_drag_asymmetry_tmp(pcols,ndir_asymmetry,begchunk:endchunk) ) + allocate( oro_drag_efflength_tmp(pcols,ndir_efflength,begchunk:endchunk) ) + oro_drag_convexity_tmp(:,:) = nan + oro_drag_asymmetry_tmp(:,:,:) = nan + oro_drag_efflength_tmp(:,:,:) = nan + + ! Read special orographic shape fields from topo file + call cam_grid_get_dim_names(grid_id, dim1name, dim2name) + call oro_drag_open_topo_file() + + found=.false. + call infld( 'OC', topo_file_ncid, dim1name, dim2name, 1, pcols, & + begchunk, endchunk, oro_drag_convexity_tmp(:,:), found, gridname='physgrid') + if(.not. found) call endrun('ERROR - oro_drag_init: topo file read error - OC') + + found=.false. + call infld( 'OA', topo_file_ncid, dim1name, 'ndir_asymmetry', dim2name, 1, pcols, 1, ndir_asymmetry, & + begchunk, endchunk, oro_drag_asymmetry_tmp(:,:,:), found, gridname='physgrid') + if(.not. found) call endrun('ERROR - oro_drag_init: topo file read error - OA') + + found=.false. + call infld( 'OL', topo_file_ncid, dim1name, 'ndir_efflength', dim2name, 1, pcols, 1, ndir_efflength, & + begchunk, endchunk, oro_drag_efflength_tmp(:,:,:), found, gridname='physgrid') + if(.not. found) call endrun('ERROR - oro_drag_init: topo file read error - OL') + + call oro_drag_close_topo_file() + + ! copy the oro drag data in pbuf + do c=begchunk,endchunk + pbuf_chunk => pbuf_get_chunk(pbuf2d, c) + call pbuf_set_field(pbuf_chunk, oro_drag_convexity_idx, oro_drag_convexity_tmp(:,c) ) + call pbuf_set_field(pbuf_chunk, oro_drag_asymmetry_idx, oro_drag_asymmetry_tmp(:,:,c) ) + call pbuf_set_field(pbuf_chunk, oro_drag_efflength_idx, oro_drag_efflength_tmp(:,:,c) ) + end do + + deallocate(oro_drag_convexity_tmp) + deallocate(oro_drag_asymmetry_tmp) + deallocate(oro_drag_efflength_tmp) + +end subroutine oro_drag_init +!========================================================================== + subroutine oro_drag_interface(state, cam_in, sgh, pbuf, dtime, nm, & gwd_ls, gwd_bl, gwd_ss, gwd_fd, & od_ls_ncleff, od_bl_ncd,od_ss_sncleff, & @@ -33,13 +169,12 @@ subroutine oro_drag_interface(state, cam_in, sgh, pbuf, dtime, nm, dusfc_ls, dvsfc_ls ,dusfc_bl, dvsfc_bl, & dusfc_ss, dvsfc_ss ,dusfc_fd, dvsfc_fd) use physics_types, only: physics_state - use physics_buffer, only: physics_buffer_desc, pbuf_get_field, pbuf_get_index use camsrfexch, only: cam_in_t use ppgrid, only: pcols,pver,pverp use physconst, only: gravit,rair,cpair,rh2o,zvir,pi use hycoef, only: etamid - - type(physics_state), intent(in) :: state ! physics state structure ! Standard deviation of orography. + !----------------------------------------------------------------------- + type(physics_state), intent(in) :: state ! physics state structure type(cam_in_t), intent(in) :: cam_in real(r8), intent(in) :: sgh(pcols) type(physics_buffer_desc), pointer :: pbuf(:) ! Physics buffer @@ -81,17 +216,21 @@ subroutine oro_drag_interface(state, cam_in, sgh, pbuf, dtime, nm, real(r8) :: zmid(pcols,pver) ! middle interface height asl (m) real(r8) :: dz(pcols,pver) ! model layer height - !real(r8) :: g - !pblh input integer :: pblh_idx = 0 integer :: kpbl2d_in(pcols) integer :: kpbl2d_reverse_in(pcols) real(r8), pointer :: pblh(:) real(r8) :: dx(pcols),dy(pcols) - !needed index + + real(r8), pointer :: oro_drag_convexity(:) + real(r8), pointer :: oro_drag_asymmetry(:,:) + real(r8), pointer :: oro_drag_efflength(:,:) + real(r8), pointer :: oro_drag_ribulk(:) ! pbuf pointer for bulk richardson number + integer :: ncol integer :: i integer :: k + !----------------------------------------------------------------------- ncol=state%ncol !convert heights above surface to heights above sea level @@ -120,7 +259,7 @@ subroutine oro_drag_interface(state, cam_in, sgh, pbuf, dtime, nm, do k=1,pver do i=1,ncol - ! assign values for level top/bottom + ! assign values for level top/bottom ztop(i,k)=state%zi(i,k) zbot(i,k)=state%zi(i,k+1) enddo @@ -145,34 +284,41 @@ subroutine oro_drag_interface(state, cam_in, sgh, pbuf, dtime, nm, kpbl2d_reverse_in(i)=pverp-kpbl2d_in(i)!pverp-k end do + call pbuf_get_field(pbuf, oro_drag_convexity_idx, oro_drag_convexity ) + call pbuf_get_field(pbuf, oro_drag_asymmetry_idx, oro_drag_asymmetry ) + call pbuf_get_field(pbuf, oro_drag_efflength_idx, oro_drag_efflength ) + call pbuf_get_field(pbuf, oro_drag_ribulk_idx, oro_drag_ribulk) + !get grid size for dx,dy call grid_size(state,dx,dy) + !interface for orographic drag - call od_gsd(& - u3d=state%u(:ncol,pver:1:-1),v3d=state%v(:ncol,pver:1:-1),t3d=state%t(:ncol,pver:1:-1),& - qv3d=state%q(:ncol,pver:1:-1,1),p3d=state%pmid(:ncol,pver:1:-1),p3di=state%pint(:ncol,pver+1:1:-1),& - pi3d=state%exner(:ncol,pver:1:-1),z=zbot(:ncol,pver:1:-1),& - od_ls_ncleff=od_ls_ncleff,od_bl_ncd=od_bl_ncd,od_ss_sncleff=od_ss_sncleff,& - rublten=utgw(:ncol,pver:1:-1),rvblten=vtgw(:ncol,pver:1:-1),rthblten=ttgw(:ncol,pver:1:-1),& - dtaux3d_ls=dtaux3_ls(:ncol,pver:1:-1),dtauy3d_ls=dtauy3_ls(:ncol,pver:1:-1),& - dtaux3d_bl=dtaux3_bl(:ncol,pver:1:-1),dtauy3d_bl=dtauy3_bl(:ncol,pver:1:-1),& - dtaux3d_ss=dtaux3_ss(:ncol,pver:1:-1),dtauy3d_ss=dtauy3_ss(:ncol,pver:1:-1),& - dtaux3d_fd=dtaux3_fd(:ncol,pver:1:-1),dtauy3d_fd=dtauy3_fd(:ncol,pver:1:-1),& - dusfcg_ls=dusfc_ls(:ncol),dvsfcg_ls=dvsfc_ls(:ncol),& - dusfcg_bl=dusfc_bl(:ncol),dvsfcg_bl=dvsfc_bl(:ncol),& - dusfcg_ss=dusfc_ss(:ncol),dvsfcg_ss=dvsfc_ss(:ncol),& - dusfcg_fd=dusfc_fd(:ncol),dvsfcg_fd=dvsfc_fd(:ncol),& - xland=cam_in%landfrac,br=state%ribulk(:ncol),& - var2d=sgh(:ncol),oc12d=state%oc(:ncol),& - oa2d=state%oadir(:ncol,:),ol2d=state%ol(:ncol,:),& - znu=etamid(pver:1:-1),dz=dz(:ncol,pver:1:-1),pblh=pblh(:ncol),& - cp=cpair,g=gravit,rd=rair,rv=rh2o,ep1=zvir,pi=pi,bnvbg=nm(:ncol,pver:1:-1),& - dt=dtime,dx=dx,dy=dy,& - kpbl2d=kpbl2d_reverse_in,gwd_opt=0,& - ids=1,ide=ncol,jds=0,jde=0,kds=1,kde=pver, & - ims=1,ime=ncol,jms=0,jme=0,kms=1,kme=pver, & - its=1,ite=ncol,jts=0,jte=0,kts=1,kte=pver, & - gwd_ls=gwd_ls,gwd_bl=gwd_bl,gwd_ss=gwd_ss,gwd_fd=gwd_fd ) + call od_gsd(u3d=state%u(:ncol,pver:1:-1),v3d=state%v(:ncol,pver:1:-1),t3d=state%t(:ncol,pver:1:-1),& + qv3d=state%q(:ncol,pver:1:-1,1),p3d=state%pmid(:ncol,pver:1:-1),p3di=state%pint(:ncol,pver+1:1:-1),& + pi3d=state%exner(:ncol,pver:1:-1),z=zbot(:ncol,pver:1:-1),& + od_ls_ncleff=od_ls_ncleff,od_bl_ncd=od_bl_ncd,od_ss_sncleff=od_ss_sncleff,& + rublten=utgw(:ncol,pver:1:-1),rvblten=vtgw(:ncol,pver:1:-1),rthblten=ttgw(:ncol,pver:1:-1),& + dtaux3d_ls=dtaux3_ls(:ncol,pver:1:-1),dtauy3d_ls=dtauy3_ls(:ncol,pver:1:-1),& + dtaux3d_bl=dtaux3_bl(:ncol,pver:1:-1),dtauy3d_bl=dtauy3_bl(:ncol,pver:1:-1),& + dtaux3d_ss=dtaux3_ss(:ncol,pver:1:-1),dtauy3d_ss=dtauy3_ss(:ncol,pver:1:-1),& + dtaux3d_fd=dtaux3_fd(:ncol,pver:1:-1),dtauy3d_fd=dtauy3_fd(:ncol,pver:1:-1),& + dusfcg_ls=dusfc_ls(:ncol),dvsfcg_ls=dvsfc_ls(:ncol),& + dusfcg_bl=dusfc_bl(:ncol),dvsfcg_bl=dvsfc_bl(:ncol),& + dusfcg_ss=dusfc_ss(:ncol),dvsfcg_ss=dvsfc_ss(:ncol),& + dusfcg_fd=dusfc_fd(:ncol),dvsfcg_fd=dvsfc_fd(:ncol),& + xland=cam_in%landfrac,br=oro_drag_ribulk(:ncol),& + var2d=sgh(:ncol),& + oc12d=oro_drag_convexity(:ncol),& + oa2d=oro_drag_asymmetry(:ncol,:),& + ol2d=oro_drag_efflength(:ncol,:),& + znu=etamid(pver:1:-1),dz=dz(:ncol,pver:1:-1),pblh=pblh(:ncol),& + cp=cpair,g=gravit,rd=rair,rv=rh2o,ep1=zvir,pi=pi,bnvbg=nm(:ncol,pver:1:-1),& + dt=dtime,dx=dx,dy=dy,& + kpbl2d=kpbl2d_reverse_in,gwd_opt=0,& + ids=1,ide=ncol,jds=0,jde=0,kds=1,kde=pver, & + ims=1,ime=ncol,jms=0,jme=0,kms=1,kme=pver, & + its=1,ite=ncol,jts=0,jte=0,kts=1,kte=pver, & + gwd_ls=gwd_ls,gwd_bl=gwd_bl,gwd_ss=gwd_ss,gwd_fd=gwd_fd ) end subroutine oro_drag_interface @@ -263,39 +409,39 @@ subroutine grid_size(state, grid_dx, grid_dy) integer :: i do i=1,state%ncol - ! determine the column area in radians - column_area = get_area_p(state%lchnk,i) - ! convert to degrees - degree = sqrt(column_area)*(180._r8/shr_const_pi) - - ! convert latitude to radians - lat_in_rad = state%lat(i)*(shr_const_pi/180._r8) - - ! Now find meters per degree latitude - ! Below equation finds distance between two points on an ellipsoid, derived from expansion - ! taking into account ellipsoid using World Geodetic System (WGS84) reference - mpdeglat = earth_ellipsoid1 - earth_ellipsoid2 * cos(2._r8*lat_in_rad) + earth_ellipsoid3 * cos(4._r8*lat_in_rad) - grid_dx(i) = mpdeglat * degree - grid_dy(i) = grid_dx(i) ! Assume these are the same + ! determine the column area in radians + column_area = get_area_p(state%lchnk,i) + ! convert to degrees + degree = sqrt(column_area)*(180._r8/shr_const_pi) + + ! convert latitude to radians + lat_in_rad = state%lat(i)*(shr_const_pi/180._r8) + + ! Now find meters per degree latitude + ! Below equation finds distance between two points on an ellipsoid, derived from expansion + ! taking into account ellipsoid using World Geodetic System (WGS84) reference + mpdeglat = earth_ellipsoid1 - earth_ellipsoid2 * cos(2._r8*lat_in_rad) + earth_ellipsoid3 * cos(4._r8*lat_in_rad) + grid_dx(i) = mpdeglat * degree + grid_dy(i) = grid_dx(i) ! Assume these are the same enddo end subroutine grid_size !========================================================================== subroutine od_gsd(u3d,v3d,t3d,qv3d,p3d,p3di,pi3d,z, & - od_ls_ncleff,od_bl_ncd,od_ss_sncleff, & - rublten,rvblten,rthblten, & - dtaux3d_ls,dtauy3d_ls,dtaux3d_bl,dtauy3d_bl, & - dtaux3d_ss,dtauy3d_ss,dtaux3d_fd,dtauy3d_fd, & - dusfcg_ls,dvsfcg_ls,dusfcg_bl,dvsfcg_bl,dusfcg_ss,dvsfcg_ss, & - dusfcg_fd,dvsfcg_fd,xland,br, & - var2d,oc12d,oa2d,ol2d,znu,znw,p_top,dz,pblh, & - cp,g,rd,rv,ep1,pi,bnvbg, & - dt,dx,dy,kpbl2d,gwd_opt, & - ids,ide, jds,jde, kds,kde, & - ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte, & - gwd_ls,gwd_bl,gwd_ss,gwd_fd) + od_ls_ncleff,od_bl_ncd,od_ss_sncleff, & + rublten,rvblten,rthblten, & + dtaux3d_ls,dtauy3d_ls,dtaux3d_bl,dtauy3d_bl, & + dtaux3d_ss,dtauy3d_ss,dtaux3d_fd,dtauy3d_fd, & + dusfcg_ls,dvsfcg_ls,dusfcg_bl,dvsfcg_bl,dusfcg_ss,dvsfcg_ss, & + dusfcg_fd,dvsfcg_fd,xland,br, & + var2d,oc12d,oa2d,ol2d,znu,znw,p_top,dz,pblh, & + cp,g,rd,rv,ep1,pi,bnvbg, & + dt,dx,dy,kpbl2d,gwd_opt, & + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte, & + gwd_ls,gwd_bl,gwd_ss,gwd_fd) !------------------------------------------------------------------------------- implicit none !------------------------------------------------------------------------------- @@ -375,8 +521,8 @@ subroutine od_gsd(u3d,v3d,t3d,qv3d,p3d,p3di,pi3d,z, & !input topographic parameters real(r8), dimension( ims:ime ), intent(in), optional :: var2d real(r8), dimension( ims:ime ), intent(in), optional :: oc12d - real(r8), dimension( ims:ime,nvar_dirOL ),intent(in), optional :: ol2d - real(r8), dimension( ims:ime,nvar_dirOA ),intent(in), optional :: oa2d + real(r8), dimension( ims:ime,ndir_efflength ),intent(in), optional :: ol2d + real(r8), dimension( ims:ime,ndir_asymmetry ),intent(in), optional :: oa2d !input model parameters real(r8), intent(in), optional :: dt real(r8), intent(in), optional :: p_top @@ -426,8 +572,8 @@ subroutine od_gsd(u3d,v3d,t3d,qv3d,p3d,p3di,pi3d,z, & real(r8), dimension( its:ite, kts:kte ) :: delprsi real(r8), dimension( its:ite, kts:kte ) :: pdh real(r8), dimension( its:ite, kts:kte+1 ) :: pdhi - real(r8), dimension( its:ite, nvar_dirOA ) :: oa4 - real(r8), dimension( its:ite, nvar_dirOL ) :: ol4 + real(r8), dimension( its:ite, ndir_asymmetry ) :: oa4 + real(r8), dimension( its:ite, ndir_efflength ) :: ol4 integer :: i,j,k,kpblmax !determine the lowest level for planet boundary layer do k = kts,kte @@ -454,58 +600,58 @@ subroutine od_gsd(u3d,v3d,t3d,qv3d,p3d,p3di,pi3d,z, & ol4(i,:) = ol2d(i,:) enddo endif - !call the od2d for calculatino of each grid - call od2d(dudt=rublten(ims,kms),dvdt=rvblten(ims,kms) & - ,dthdt=rthblten(ims,kms) & - ,ncleff=od_ls_ncleff,ncd=od_bl_ncd,sncleff=od_ss_sncleff & - ,dtaux2d_ls=dtaux2d_ls,dtauy2d_ls=dtauy2d_ls & - ,dtaux2d_bl=dtaux2d_bl,dtauy2d_bl=dtauy2d_bl & - ,dtaux2d_ss=dtaux2d_ss,dtauy2d_ss=dtauy2d_ss & - ,dtaux2d_fd=dtaux2d_fd,dtauy2d_fd=dtauy2d_fd & - ,u1=u3d(ims,kms),v1=v3d(ims,kms) & - ,t1=t3d(ims,kms) & - ,q1=qv3d(ims,kms) & - ,del=delprsi(its,kts) & - ,prsi=pdhi(its,kts) & - ,prsl=pdh(its,kts),prslk=pi3d(ims,kms) & - ,zl=z(ims,kms),rcl=1.0_r8 & - ,xland1=xland(ims),br1=br(ims),hpbl=pblh(ims) & - ,bnv_in=bnvbg(ims,kms) & - ,dz2=dz(ims,kms) & - ,kpblmax=kpblmax & - ,dusfc_ls=dusfc_ls,dvsfc_ls=dvsfc_ls & - ,dusfc_bl=dusfc_bl,dvsfc_bl=dvsfc_bl & - ,dusfc_ss=dusfc_ss,dvsfc_ss=dvsfc_ss & - ,dusfc_fd=dusfc_fd,dvsfc_fd=dvsfc_fd & - ,var=var2d(ims),oc1=oc12d(ims) & - ,oa4=oa4,ol4=ol4 & - ,g=g,cp=cp,rd=rd,rv=rv,fv=ep1,pi=pi & - ,dxmeter=dx,dymeter=dy,deltim=dt & - ,kpbl=kpbl2d(ims) & - ,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde & - ,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme & - ,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte & - ,gsd_gwd_ls=gwd_ls,gsd_gwd_bl=gwd_bl,gsd_gwd_ss=gwd_ss,gsd_gwd_fd=gwd_fd) - !set the total stress output to each terms for the 4 drag schemes - do i = its,ite - dusfcg_ls(i)=dusfc_ls(i) - dvsfcg_ls(i)=dvsfc_ls(i) - dusfcg_bl(i)=dusfc_bl(i) - dvsfcg_bl(i)=dvsfc_bl(i) - dusfcg_ss(i)=dusfc_ss(i) - dvsfcg_ss(i)=dvsfc_ss(i) - dusfcg_fd(i)=dusfc_fd(i) - dvsfcg_fd(i)=dvsfc_fd(i) - enddo - !set the 3D output tendencies to each terms for the 4 drag schemes - dtaux3d_ls=dtaux2d_ls - dtaux3d_bl=dtaux2d_bl - dtauy3d_ls=dtauy2d_ls - dtauy3d_bl=dtauy2d_bl - dtaux3d_ss=dtaux2d_ss - dtaux3d_fd=dtaux2d_fd - dtauy3d_ss=dtauy2d_ss - dtauy3d_fd=dtauy2d_fd + !call the od2d for calculatino of each grid + call od2d(dudt=rublten(ims,kms),dvdt=rvblten(ims,kms) & + ,dthdt=rthblten(ims,kms) & + ,ncleff=od_ls_ncleff,ncd=od_bl_ncd,sncleff=od_ss_sncleff & + ,dtaux2d_ls=dtaux2d_ls,dtauy2d_ls=dtauy2d_ls & + ,dtaux2d_bl=dtaux2d_bl,dtauy2d_bl=dtauy2d_bl & + ,dtaux2d_ss=dtaux2d_ss,dtauy2d_ss=dtauy2d_ss & + ,dtaux2d_fd=dtaux2d_fd,dtauy2d_fd=dtauy2d_fd & + ,u1=u3d(ims,kms),v1=v3d(ims,kms) & + ,t1=t3d(ims,kms) & + ,q1=qv3d(ims,kms) & + ,del=delprsi(its,kts) & + ,prsi=pdhi(its,kts) & + ,prsl=pdh(its,kts),prslk=pi3d(ims,kms) & + ,zl=z(ims,kms),rcl=1.0_r8 & + ,xland1=xland(ims),br1=br(ims),hpbl=pblh(ims) & + ,bnv_in=bnvbg(ims,kms) & + ,dz2=dz(ims,kms) & + ,kpblmax=kpblmax & + ,dusfc_ls=dusfc_ls,dvsfc_ls=dvsfc_ls & + ,dusfc_bl=dusfc_bl,dvsfc_bl=dvsfc_bl & + ,dusfc_ss=dusfc_ss,dvsfc_ss=dvsfc_ss & + ,dusfc_fd=dusfc_fd,dvsfc_fd=dvsfc_fd & + ,var=var2d(ims),oc1=oc12d(ims) & + ,oa4=oa4,ol4=ol4 & + ,g=g,cp=cp,rd=rd,rv=rv,fv=ep1,pi=pi & + ,dxmeter=dx,dymeter=dy,deltim=dt & + ,kpbl=kpbl2d(ims) & + ,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde & + ,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme & + ,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte & + ,gsd_gwd_ls=gwd_ls,gsd_gwd_bl=gwd_bl,gsd_gwd_ss=gwd_ss,gsd_gwd_fd=gwd_fd) + !set the total stress output to each terms for the 4 drag schemes + do i = its,ite + dusfcg_ls(i) = dusfc_ls(i) + dvsfcg_ls(i) = dvsfc_ls(i) + dusfcg_bl(i) = dusfc_bl(i) + dvsfcg_bl(i) = dvsfc_bl(i) + dusfcg_ss(i) = dusfc_ss(i) + dvsfcg_ss(i) = dvsfc_ss(i) + dusfcg_fd(i) = dusfc_fd(i) + dvsfcg_fd(i) = dvsfc_fd(i) + enddo + !set the 3D output tendencies to each terms for the 4 drag schemes + dtaux3d_ls = dtaux2d_ls + dtaux3d_bl = dtaux2d_bl + dtauy3d_ls = dtauy2d_ls + dtauy3d_bl = dtauy2d_bl + dtaux3d_ss = dtaux2d_ss + dtaux3d_fd = dtaux2d_fd + dtauy3d_ss = dtauy2d_ss + dtauy3d_fd = dtauy2d_fd end subroutine od_gsd ! @@ -513,21 +659,21 @@ end subroutine od_gsd ! !------------------------------------------------------------------------------- subroutine od2d(dudt,dvdt,dthdt,ncleff,ncd,sncleff, & - dtaux2d_ls,dtauy2d_ls, & - dtaux2d_bl,dtauy2d_bl, & - dtaux2d_ss,dtauy2d_ss, & - dtaux2d_fd,dtauy2d_fd, & - u1,v1,t1,q1, & - del, & - prsi,prsl,prslk,zl,rcl, & - xland1,br1,hpbl,bnv_in,dz2, & - kpblmax,dusfc_ls,dvsfc_ls,dusfc_bl,dvsfc_bl, & - dusfc_ss,dvsfc_ss,dusfc_fd,dvsfc_fd,var,oc1,oa4,ol4, & - g,cp,rd,rv,fv,pi,dxmeter,dymeter,deltim,kpbl, & - ids,ide, jds,jde, kds,kde, & - ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte, & - gsd_gwd_ls,gsd_gwd_bl,gsd_gwd_ss,gsd_gwd_fd) + dtaux2d_ls,dtauy2d_ls, & + dtaux2d_bl,dtauy2d_bl, & + dtaux2d_ss,dtauy2d_ss, & + dtaux2d_fd,dtauy2d_fd, & + u1,v1,t1,q1, & + del, & + prsi,prsl,prslk,zl,rcl, & + xland1,br1,hpbl,bnv_in,dz2, & + kpblmax,dusfc_ls,dvsfc_ls,dusfc_bl,dvsfc_bl, & + dusfc_ss,dvsfc_ss,dusfc_fd,dvsfc_fd,var,oc1,oa4,ol4, & + g,cp,rd,rv,fv,pi,dxmeter,dymeter,deltim,kpbl, & + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte, & + gsd_gwd_ls,gsd_gwd_bl,gsd_gwd_ss,gsd_gwd_fd) ! This code handles the time tendencies of u v due to the effect of mountain ! induced gravity wave drag from sub-grid scale orography. It includes 4 parts: ! orographic gravity wave drag and flow-blocking drag (Xie et al.,2020),small-scale @@ -582,8 +728,8 @@ subroutine od2d(dudt,dvdt,dthdt,ncleff,ncd,sncleff, & real(r8), dimension(:), intent(in) :: dxmeter real(r8), dimension(:), intent(in) :: dymeter !input topo variables - real(r8), dimension( ims:ime,nvar_dirOA ), intent(in) :: oa4 - real(r8), dimension( ims:ime,nvar_dirOL ), intent(in) :: ol4 + real(r8), dimension( ims:ime,ndir_asymmetry ), intent(in) :: oa4 + real(r8), dimension( ims:ime,ndir_efflength ), intent(in) :: ol4 real(r8), dimension( ims:ime ) , intent(in) :: var real(r8), dimension( ims:ime ) , intent(in) :: oc1 !input atmospheric variables @@ -673,7 +819,7 @@ subroutine od2d(dudt,dvdt,dthdt,ncleff,ncd,sncleff, & real(r8),parameter :: tndmax = 400._r8 / 86400._r8 ! convert 400 m/s/day to m/s/s integer,parameter :: kpblmin = 2 !number of direction for ogwd - integer,parameter :: mdir=2*nvar_dirOL + integer,parameter :: mdir=2*ndir_efflength ! variables for flow-blocking drag real(r8),parameter :: frmax = 10._r8 real(r8),parameter :: olmin = 1.0e-5_r8 @@ -736,8 +882,8 @@ subroutine od2d(dudt,dvdt,dthdt,ncleff,ncd,sncleff, & real(r8),dimension( its:ite ) :: dely real(r8),dimension( its:ite ) :: dxy real(r8),dimension( its:ite ) :: dxyp - real(r8),dimension( its:ite,nvar_dirOL ):: dxy4 - real(r8),dimension( its:ite,nvar_dirOL ):: dxy4p + real(r8),dimension( its:ite,ndir_efflength ):: dxy4 + real(r8),dimension( its:ite,ndir_efflength ):: dxy4p !topo parameters real(r8),dimension( its:ite ) :: olp real(r8),dimension( its:ite ) :: od @@ -893,192 +1039,192 @@ subroutine od2d(dudt,dvdt,dthdt,ncleff,ncd,sncleff, & ! ! For ls and bl only IF (gsd_gwd_ls.or.gsd_gwd_bl) then - ! figure out low-level horizontal wind direction - ! order into a counterclockwise index instead - ! - do i = its,ite - wdir = atan2(vbar(i),ubar(i)) + pi!changed into y/x - wdir1 = wdir-pi - if (wdir1.ge.0._r8.and.wdir1.lt.pi) then - nwd = MOD(nint(fdir*wdir1),mdir) + 1 - else!(-pi,0) - nwd = MOD(nint(fdir*(wdir1+2._r8*pi)),mdir) + 1 - endif - !turn backwords because start is pi - !need turning - rad = 4.0_r8*atan(1.0_r8)/180.0_r8 - theta = (real(nwd,kind=r8)-1._r8)*(360._r8/real(mdir,kind=r8)) - !select OA - oa1(i) = oa4(i,1)*cos(theta*rad)+oa4(i,2)*sin(theta*rad) - !select OL - ol(i) = ol4(i,MOD(nwd-1,int(mdir/2))+1) - !calculate dxygrid, not so slow - call dxygrid(dxmeter(i),dymeter(i),theta,dxy(i)) + ! figure out low-level horizontal wind direction + ! order into a counterclockwise index instead ! - !----- compute orographic width along (ol) and perpendicular (olp) - !----- the direction of wind - !put wdir inside the (0,2*pi) section - !changing pi/2 either way is perpendicular - !wdir1=wdir-pi + do i = its,ite + wdir = atan2(vbar(i),ubar(i)) + pi!changed into y/x + wdir1 = wdir-pi + if (wdir1.ge.0._r8.and.wdir1.lt.pi) then + nwd = MOD(nint(fdir*wdir1),mdir) + 1 + else!(-pi,0) + nwd = MOD(nint(fdir*(wdir1+2._r8*pi)),mdir) + 1 + endif + !turn backwords because start is pi + !need turning + rad = 4.0_r8*atan(1.0_r8)/180.0_r8 + theta = (real(nwd,kind=r8)-1._r8)*(360._r8/real(mdir,kind=r8)) + !select OA + oa1(i) = oa4(i,1)*cos(theta*rad)+oa4(i,2)*sin(theta*rad) + !select OL + ol(i) = ol4(i,MOD(nwd-1,int(mdir/2))+1) + !calculate dxygrid, not so slow + call dxygrid(dxmeter(i),dymeter(i),theta,dxy(i)) + ! + !----- compute orographic width along (ol) and perpendicular (olp) + !----- the direction of wind + !put wdir inside the (0,2*pi) section + !changing pi/2 either way is perpendicular + !wdir1=wdir-pi if (wdir1.ge.0._r8.and.wdir1.lt.pi) then - nwd1 = MOD(nint(fdir*(wdir1+pi/2._r8)),mdir) + 1 - olp(i)=ol4(i,MOD(nwd1-1,int(mdir/2))+1) + nwd1 = MOD(nint(fdir*(wdir1+pi/2._r8)),mdir) + 1 + olp(i)=ol4(i,MOD(nwd1-1,int(mdir/2))+1) else!(-pi,0) - nwd1 = MOD(nint(fdir*(wdir1-pi/2._r8+2._r8*pi)),mdir) + 1 - olp(i)=ol4(i,MOD(nwd1-1,int(mdir/2))+1) + nwd1 = MOD(nint(fdir*(wdir1-pi/2._r8+2._r8*pi)),mdir) + 1 + olp(i)=ol4(i,MOD(nwd1-1,int(mdir/2))+1) endif theta=(real(nwd1,kind=r8)-1._r8)*(360._r8/real(mdir,kind=r8)) call dxygrid(dxmeter(i),dymeter(i),theta,dxyp(i)) - ! - ! - !----- compute orographic direction (horizontal orographic aspect ratio) - ! - od(i) = olp(i)/max(ol(i),olmin) - od(i) = min(od(i),odmax) - od(i) = max(od(i),odmin) - ! - !----- compute length of grid in the along(dxy) and cross(dxyp) wind directions - ! - enddo + ! + ! + !----- compute orographic direction (horizontal orographic aspect ratio) + ! + od(i) = olp(i)/max(ol(i),olmin) + od(i) = min(od(i),odmax) + od(i) = max(od(i),odmin) + ! + !----- compute length of grid in the along(dxy) and cross(dxyp) wind directions + ! + enddo ENDIF !============================================ ! END INITIALIZATION; BEGIN GWD CALCULATIONS: !============================================ IF (gsd_gwd_ls.or.gsd_gwd_bl.and.(ls_taper .GT. 1.E-02) ) THEN - ! - !--- saving richardson number in usqj for migwdi - ! - do k = kts,kte-1 + ! + !--- saving richardson number in usqj for migwdi + ! + do k = kts,kte-1 + do i = its,ite + ti = 2.0_r8 / (t1(i,k)+t1(i,k+1)) + rdz = 1._r8/(zl(i,k+1) - zl(i,k)) + tem1 = u1(i,k) - u1(i,k+1) + tem2 = v1(i,k) - v1(i,k+1) + dw2 = rcl*(tem1*tem1 + tem2*tem2) + shr2 = max(dw2,dw2min) * rdz * rdz + bvf2 = g*(g/cp+rdz*(vtj(i,k+1)-vtj(i,k))) * ti + usqj(i,k) = max(bvf2/shr2,rimin) + bnv2(i,k) = max(bnv_in(i,k)**2,bnv2min ) + enddo + enddo + ! + !----compute the "low level" or 1/3 wind magnitude (m/s) + ! do i = its,ite - ti = 2.0_r8 / (t1(i,k)+t1(i,k+1)) - rdz = 1._r8/(zl(i,k+1) - zl(i,k)) - tem1 = u1(i,k) - u1(i,k+1) - tem2 = v1(i,k) - v1(i,k+1) - dw2 = rcl*(tem1*tem1 + tem2*tem2) - shr2 = max(dw2,dw2min) * rdz * rdz - bvf2 = g*(g/cp+rdz*(vtj(i,k+1)-vtj(i,k))) * ti - usqj(i,k) = max(bvf2/shr2,rimin) - bnv2(i,k) = max(bnv_in(i,k)**2,bnv2min ) + ulow(i) = max(sqrt(ubar(i)*ubar(i) + vbar(i)*vbar(i)), 1.0_r8) + rulow(i) = 1._r8/ulow(i) enddo - enddo - ! - !----compute the "low level" or 1/3 wind magnitude (m/s) - ! - do i = its,ite - ulow(i) = max(sqrt(ubar(i)*ubar(i) + vbar(i)*vbar(i)), 1.0_r8) - rulow(i) = 1._r8/ulow(i) - enddo - do k = kts,kte-1 + do k = kts,kte-1 + do i = its,ite + velco(i,k) = (0.5_r8*rcs) * ((u1(i,k)+u1(i,k+1)) * ubar(i) & + + (v1(i,k)+v1(i,k+1)) * vbar(i)) + velco(i,k) = velco(i,k) * rulow(i) + if ((velco(i,k).lt.veleps) .and. (velco(i,k).gt.0._r8)) then + velco(i,k) = veleps + endif + enddo + enddo + ! + ! no drag when critical level in the base layer + ! do i = its,ite - velco(i,k) = (0.5_r8*rcs) * ((u1(i,k)+u1(i,k+1)) * ubar(i) & - + (v1(i,k)+v1(i,k+1)) * vbar(i)) - velco(i,k) = velco(i,k) * rulow(i) - if ((velco(i,k).lt.veleps) .and. (velco(i,k).gt.0._r8)) then - velco(i,k) = veleps - endif + ldrag(i) = velco(i,1).le.0._r8 enddo - enddo - ! - ! no drag when critical level in the base layer - ! - do i = its,ite - ldrag(i) = velco(i,1).le.0._r8 - enddo - ! - ! no drag when velco.lt.0 - ! - do k = kpblmin,kpblmax + ! + ! no drag when velco.lt.0 + ! + do k = kpblmin,kpblmax + do i = its,ite + if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. velco(i,k).le.0._r8 + enddo + enddo + ! + ! no drag when bnv2.lt.0 + ! + do k = kts,kpblmax + do i = its,ite + if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. bnv2(i,k).lt.0._r8 + enddo + enddo + ! + !-----the low level weighted average ri is stored in usqj(1,1; im) + !-----the low level weighted average n**2 is stored in bnv2(1,1; im) + !---- this is called bnvl2 in phys_gwd_alpert_sub not bnv2 + !---- rdelks (del(k)/delks) vert ave factor so we can * instead of / + ! do i = its,ite - if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. velco(i,k).le.0._r8 + wtkbj = (prsl(i,1)-prsl(i,2)) * delks1(i) + bnv2(i,1) = wtkbj * bnv2(i,1) + usqj(i,1) = wtkbj * usqj(i,1) enddo - enddo - ! - ! no drag when bnv2.lt.0 - ! - do k = kts,kpblmax + + do k = kpblmin,kpblmax + do i = its,ite + if (k .lt. kbl(i)) then + rdelks = (prsl(i,k)-prsl(i,k+1)) * delks1(i) + bnv2(i,1) = bnv2(i,1) + bnv2(i,k) * rdelks + usqj(i,1) = usqj(i,1) + usqj(i,k) * rdelks + endif + enddo + enddo + do i = its,ite - if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. bnv2(i,k).lt.0._r8 + ldrag(i) = ldrag(i) .or. bnv2(i,1).le.0.0_r8 + ldrag(i) = ldrag(i) .or. ulow(i) .eq.1.0_r8 + ldrag(i) = ldrag(i) .or. var(i) .le.0.0_r8 enddo - enddo - ! - !-----the low level weighted average ri is stored in usqj(1,1; im) - !-----the low level weighted average n**2 is stored in bnv2(1,1; im) - !---- this is called bnvl2 in phys_gwd_alpert_sub not bnv2 - !---- rdelks (del(k)/delks) vert ave factor so we can * instead of / - ! - do i = its,ite - wtkbj = (prsl(i,1)-prsl(i,2)) * delks1(i) - bnv2(i,1) = wtkbj * bnv2(i,1) - usqj(i,1) = wtkbj * usqj(i,1) - enddo - - do k = kpblmin,kpblmax + ! + ! set all ri low level values to the low level value + ! + do k = kpblmin,kpblmax + do i = its,ite + if (k .lt. kbl(i)) usqj(i,k) = usqj(i,1) + enddo + enddo + do i = its,ite - if (k .lt. kbl(i)) then - rdelks = (prsl(i,k)-prsl(i,k+1)) * delks1(i) - bnv2(i,1) = bnv2(i,1) + bnv2(i,k) * rdelks - usqj(i,1) = usqj(i,1) + usqj(i,k) * rdelks + if (.not.ldrag(i)) then + bnv(i) = sqrt( bnv2(i,1) ) + fr(i) = bnv(i) * rulow(i) * 2._r8 * var(i) * od(i) + fr(i) = min(fr(i),frmax) + xn(i) = ubar(i) * rulow(i) + yn(i) = vbar(i) * rulow(i) endif enddo - enddo - - do i = its,ite - ldrag(i) = ldrag(i) .or. bnv2(i,1).le.0.0_r8 - ldrag(i) = ldrag(i) .or. ulow(i) .eq.1.0_r8 - ldrag(i) = ldrag(i) .or. var(i) .le.0.0_r8 - enddo - ! - ! set all ri low level values to the low level value - ! - do k = kpblmin,kpblmax + ! + ! compute the base level stress and store it in taub + ! calculate enhancement factor, number of mountains & aspect + ! ratio const. use simplified relationship between standard + ! deviation & critical hgt + ! do i = its,ite - if (k .lt. kbl(i)) usqj(i,k) = usqj(i,1) + if (.not. ldrag(i)) then + !maintain (oa+2) greater than or equal to 0 + efact = max(oa1(i)+2._r8,0._r8) ** (ce*fr(i)/frc) + efact = min(max(efact,efmin),efmax) + ! cleff (effective grid length) is highly tunable parameter + ! the bigger (smaller) value produce weaker (stronger) wave drag + cleff = sqrt(dxy(i)**2._r8 + dxyp(i)**2._r8) + !tune the times of drag + cleff = (3._r8/ncleff) * max(dxmax_ls,cleff) + coefm(i) = (1._r8 + ol(i)) ** (oa1(i)+1._r8) + xlinv(i) = coefm(i) / cleff + tem = fr(i) * fr(i) * 1.!oc1(i) + gfobnv = gmax * tem / ((tem + cg)*bnv(i)) + ! + if (gsd_gwd_ls) then + taub(i) = xlinv(i) * roll(i) * ulow(i) * ulow(i) & + * ulow(i) * gfobnv * efact + else ! We've gotten what we need for the blocking scheme + taub(i) = 0.0_r8 + end if + else + taub(i) = 0.0_r8 + xn(i) = 0.0_r8 + yn(i) = 0.0_r8 + endif enddo - enddo - - do i = its,ite - if (.not.ldrag(i)) then - bnv(i) = sqrt( bnv2(i,1) ) - fr(i) = bnv(i) * rulow(i) * 2._r8 * var(i) * od(i) - fr(i) = min(fr(i),frmax) - xn(i) = ubar(i) * rulow(i) - yn(i) = vbar(i) * rulow(i) - endif - enddo - ! - ! compute the base level stress and store it in taub - ! calculate enhancement factor, number of mountains & aspect - ! ratio const. use simplified relationship between standard - ! deviation & critical hgt - ! - do i = its,ite - if (.not. ldrag(i)) then - !maintain (oa+2) greater than or equal to 0 - efact = max(oa1(i)+2._r8,0._r8) ** (ce*fr(i)/frc) - efact = min(max(efact,efmin),efmax) - ! cleff (effective grid length) is highly tunable parameter - ! the bigger (smaller) value produce weaker (stronger) wave drag - cleff = sqrt(dxy(i)**2._r8 + dxyp(i)**2._r8) - !tune the times of drag - cleff = (3._r8/ncleff) * max(dxmax_ls,cleff) - coefm(i) = (1._r8 + ol(i)) ** (oa1(i)+1._r8) - xlinv(i) = coefm(i) / cleff - tem = fr(i) * fr(i) * 1.!oc1(i) - gfobnv = gmax * tem / ((tem + cg)*bnv(i)) - ! - if (gsd_gwd_ls) then - taub(i) = xlinv(i) * roll(i) * ulow(i) * ulow(i) & - * ulow(i) * gfobnv * efact - else ! We've gotten what we need for the blocking scheme - taub(i) = 0.0_r8 - end if - else - taub(i) = 0.0_r8 - xn(i) = 0.0_r8 - yn(i) = 0.0_r8 - endif - enddo ENDIF ! (gsd_gwd_ls .eq. .true.).or.(gsd_gwd_bl .eq..true.) !========================================================= @@ -1092,93 +1238,93 @@ subroutine od2d(dudt,dvdt,dthdt,ncleff,ncd,sncleff, & zq=0._r8 IF (gsd_gwd_ss.and.(ss_taper.GT.1.E-02)) THEN - ! - ! declaring potential temperature - ! - do k = kts,kte - do i = its,ite - thx(i,k) = t1(i,k)/prslk(i,k) + ! + ! declaring potential temperature + ! + do k = kts,kte + do i = its,ite + thx(i,k) = t1(i,k)/prslk(i,k) + enddo enddo - enddo - do k = kts,kte - do i = its,ite - tvcon = (1._r8+fv*q1(i,k)) - thvx(i,k) = thx(i,k)*tvcon + do k = kts,kte + do i = its,ite + tvcon = (1._r8+fv*q1(i,k)) + thvx(i,k) = thx(i,k)*tvcon + enddo enddo - enddo - ! - ! Defining layer height - ! - do k = kts,kte - do i = its,ite - zq(i,k+1) = dz2(i,k)+zq(i,k) + ! + ! Defining layer height + ! + do k = kts,kte + do i = its,ite + zq(i,k+1) = dz2(i,k)+zq(i,k) + enddo enddo - enddo - do k = kts,kte - do i = its,ite - za(i,k) = 0.5_r8*(zq(i,k)+zq(i,k+1)) + do k = kts,kte + do i = its,ite + za(i,k) = 0.5_r8*(zq(i,k)+zq(i,k+1)) + enddo enddo - enddo - do i=its,ite - hpbl2 = hpbl(i)+10._r8 - kpbl2 = kpbl(i) - kvar = 1 - do k=kts+1,MAX(kpbl(i),kts+1) - IF (za(i,k)>300._r8) then - kpbl2 = k - IF (k == kpbl(i)) then - hpbl2 = hpbl(i)+10._r8 - ELSE - hpbl2 = za(i,k)+10._r8 + do i=its,ite + hpbl2 = hpbl(i)+10._r8 + kpbl2 = kpbl(i) + kvar = 1 + do k=kts+1,MAX(kpbl(i),kts+1) + IF (za(i,k)>300._r8) then + kpbl2 = k + IF (k == kpbl(i)) then + hpbl2 = hpbl(i)+10._r8 + ELSE + hpbl2 = za(i,k)+10._r8 + ENDIF + exit ENDIF - exit - ENDIF - enddo + enddo - if(xland1(i).gt.0._r8 .and. 2._r8*var(i).le.hpbl(i))then - if(br1(i).gt.0._r8 .and. thvx(i,kpbl2)-thvx(i,kts) > 0._r8)then - cleff = sqrt(dxy(i)**2_r8 + dxyp(i)**2_r8) - cleff = (2.0_r8/sncleff) * max(dxmax_ss,cleff) - coefm(i) = (1._r8 + ol(i)) ** (oa1(i)+1._r8) - xlinv(i) = coefm(i) / cleff - govrth(i)=g/(0.5_r8*(thvx(i,kpbl2)+thvx(i,kts))) - bnrf=sqrt(govrth(i)*(thvx(i,kpbl2)-thvx(i,kts))/hpbl2) - - if(abs(bnrf/u1(i,kpbl2)).gt.xlinv(i))then - tauwavex0=0.5_r8*bnrf*xlinv(i)*(2._r8*MIN(var(i),varmax))**2_r8*ro(i,kvar)*u1(i,kvar) - tauwavex0=tauwavex0*ss_taper ! "Scale-awareness" - else - tauwavex0=0._r8 - endif + if(xland1(i).gt.0._r8 .and. 2._r8*var(i).le.hpbl(i))then + if(br1(i).gt.0._r8 .and. thvx(i,kpbl2)-thvx(i,kts) > 0._r8)then + cleff = sqrt(dxy(i)**2_r8 + dxyp(i)**2_r8) + cleff = (2.0_r8/sncleff) * max(dxmax_ss,cleff) + coefm(i) = (1._r8 + ol(i)) ** (oa1(i)+1._r8) + xlinv(i) = coefm(i) / cleff + govrth(i)=g/(0.5_r8*(thvx(i,kpbl2)+thvx(i,kts))) + bnrf=sqrt(govrth(i)*(thvx(i,kpbl2)-thvx(i,kts))/hpbl2) + + if(abs(bnrf/u1(i,kpbl2)).gt.xlinv(i))then + tauwavex0=0.5_r8*bnrf*xlinv(i)*(2._r8*MIN(var(i),varmax))**2_r8*ro(i,kvar)*u1(i,kvar) + tauwavex0=tauwavex0*ss_taper ! "Scale-awareness" + else + tauwavex0=0._r8 + endif - if(abs(bnrf/v1(i,kpbl2)).gt.xlinv(i))then - tauwavey0=0.5_r8*bnrf*xlinv(i)*(2._r8*MIN(var(i),varmax))**2._r8*ro(i,kvar)*v1(i,kvar) - tauwavey0=tauwavey0*ss_taper ! "Scale-awareness" - else - tauwavey0=0._r8 - endif + if(abs(bnrf/v1(i,kpbl2)).gt.xlinv(i))then + tauwavey0=0.5_r8*bnrf*xlinv(i)*(2._r8*MIN(var(i),varmax))**2._r8*ro(i,kvar)*v1(i,kvar) + tauwavey0=tauwavey0*ss_taper ! "Scale-awareness" + else + tauwavey0=0._r8 + endif - do k=kts,kpbl(i) !MIN(kpbl2+1,kte-1) - utendwave(i,k)=-1._r8*tauwavex0*2._r8*max((1._r8-za(i,k)/hpbl2),0._r8)/hpbl2 - vtendwave(i,k)=-1._r8*tauwavey0*2._r8*max((1._r8-za(i,k)/hpbl2),0._r8)/hpbl2 - enddo + do k=kts,kpbl(i) !MIN(kpbl2+1,kte-1) + utendwave(i,k)=-1._r8*tauwavex0*2._r8*max((1._r8-za(i,k)/hpbl2),0._r8)/hpbl2 + vtendwave(i,k)=-1._r8*tauwavey0*2._r8*max((1._r8-za(i,k)/hpbl2),0._r8)/hpbl2 + enddo + endif endif - endif - enddo ! end i loop + enddo ! end i loop - do k = kts,kte - do i = its,ite - dudt(i,k) = dudt(i,k) + utendwave(i,k) - dvdt(i,k) = dvdt(i,k) + vtendwave(i,k) - dtaux2d_ss(i,k) = utendwave(i,k) - dtauy2d_ss(i,k) = vtendwave(i,k) - dusfc_ss(i) = dusfc_ss(i) + utendwave(i,k) * del(i,k) - dvsfc_ss(i) = dvsfc_ss(i) + vtendwave(i,k) * del(i,k) + do k = kts,kte + do i = its,ite + dudt(i,k) = dudt(i,k) + utendwave(i,k) + dvdt(i,k) = dvdt(i,k) + vtendwave(i,k) + dtaux2d_ss(i,k) = utendwave(i,k) + dtauy2d_ss(i,k) = vtendwave(i,k) + dusfc_ss(i) = dusfc_ss(i) + utendwave(i,k) * del(i,k) + dvsfc_ss(i) = dvsfc_ss(i) + vtendwave(i,k) * del(i,k) + enddo enddo - enddo ENDIF ! end if gsd_gwd_ss == .true. !================================================================ @@ -1240,108 +1386,108 @@ subroutine od2d(dudt,dvdt,dthdt,ncleff,ncd,sncleff, & dvsfc_fd(i) = dvsfc_fd(i) + vtendform(i,k) * del(i,k) enddo enddo - ENDIF ! end if gsd_gwd_fd == .true. + ENDIF ! end if gsd_gwd_fd == .true. !======================================================= ! More for the large-scale gwd component !======================================================= IF (gsd_gwd_ls.and.(ls_taper.GT.1.E-02) ) THEN - ! - ! now compute vertical structure of the stress. - ! - do k = kts,kpblmax - do i = its,ite - if (k .le. kbl(i)) taup(i,k) = taub(i) - enddo - enddo - - if (scorer_on) then - ! - !determination of the interface height for scorer adjustment - ! - do i=its,ite - iint=.false. - do k=kpblmin,kte-1 - if (k.gt.kbl(i).and.usqj(i,k)-usqj(i,k-1).lt.0.and.(.not.iint)) then - iint=.true. - zl_hint(i)=zl(i,k+1) - endif + ! + ! now compute vertical structure of the stress. + ! + do k = kts,kpblmax + do i = its,ite + if (k .le. kbl(i)) taup(i,k) = taub(i) enddo enddo - endif - do k = kpblmin, kte-1 ! vertical level k loop! - kp1 = k + 1 - do i = its,ite - ! - ! unstablelayer if ri < ric - ! unstable layer if upper air vel comp along surf vel <=0 (crit lay) - ! at (u-c)=0. crit layer exists and bit vector should be set (.le.) - ! - if (k .ge. kbl(i)) then - !we modify the criteria for unstable layer - !that the lv is critical under 0.25 - !while we keep wave breaking ric for - !other larger lv - icrilv(i) = icrilv(i) .or. ( usqj(i,k) .lt. ric_rig)& - .or. (velco(i,k) .le. 0.0_r8) - brvf(i) = max(bnv2(i,k),bnv2min) ! brunt-vaisala frequency squared - brvf(i) = sqrt(brvf(i)) ! brunt-vaisala frequency - endif - enddo + if (scorer_on) then + ! + !determination of the interface height for scorer adjustment + ! + do i=its,ite + iint=.false. + do k=kpblmin,kte-1 + if (k.gt.kbl(i).and.usqj(i,k)-usqj(i,k-1).lt.0.and.(.not.iint)) then + iint=.true. + zl_hint(i)=zl(i,k+1) + endif + enddo + enddo + endif - do i = its,ite - if (k .ge. kbl(i) .and. (.not. ldrag(i))) then - if (.not.icrilv(i) .and. taup(i,k) .gt. 0.0_r8 ) then - temv = 1.0_r8 / velco(i,k) - tem1 = coefm(i)/(dxy(i)/ncleff)*(ro(i,kp1)+ro(i,k))*brvf(i)*velco(i,k)*0.5_r8 - hd = sqrt(taup(i,k) / tem1) - fro = brvf(i) * hd * temv - ! - ! rim is the minimum-richardson number by shutts (1985) - ! - tem2 = sqrt(usqj(i,k)) - tem = 1._r8 + tem2 * fro - rim = usqj(i,k) * (1._r8-fro) / (tem * tem) + do k = kpblmin, kte-1 ! vertical level k loop! + kp1 = k + 1 + do i = its,ite + ! + ! unstablelayer if ri < ric + ! unstable layer if upper air vel comp along surf vel <=0 (crit lay) + ! at (u-c)=0. crit layer exists and bit vector should be set (.le.) + ! + if (k .ge. kbl(i)) then + !we modify the criteria for unstable layer + !that the lv is critical under 0.25 + !while we keep wave breaking ric for + !other larger lv + icrilv(i) = icrilv(i) .or. ( usqj(i,k) .lt. ric_rig)& + .or. (velco(i,k) .le. 0.0_r8) + brvf(i) = max(bnv2(i,k),bnv2min) ! brunt-vaisala frequency squared + brvf(i) = sqrt(brvf(i)) ! brunt-vaisala frequency + endif + enddo - ! - ! check stability to employ the 'saturation hypothesis' - ! of lindzen (1981) except at tropospheric downstream regions - ! - if (rim .le. ric) then ! saturation hypothesis! - if ((oa1(i) .le. 0._r8).or.(kp1 .ge. kpblmin )) then - temc = 2.0_r8 + 1.0_r8 / tem2 - hd = velco(i,k) * (2.0_r8*sqrt(temc)-temc) / brvf(i) - taup(i,kp1) = tem1 * hd * hd - ! - ! taup is restricted to monotoncally decrease - ! to avoid unexpected high taup in calculation - ! - taup(i,kp1)=min(tem1*hd*hd,taup(i,k)) - ! - ! add vertical decrease at low level below hint (Kim and Doyle 2005) - ! where Ri first decreases - ! - if (scorer_on.and.k.gt.klowtop(i).and.zl(i,k).le.zl_hint(i).and.k.lt.kte-1) then - l1=(9.81_r8*bnv2(i,kp1)/velco(i,kp1)**2) - l2=(9.81_r8*bnv2(i,k)/velco(i,k)**2) - taup(i,kp1)=min(taup(i,k),taup(i,k)*(l1/l2),tem1*hd*hd) + do i = its,ite + if (k .ge. kbl(i) .and. (.not. ldrag(i))) then + if (.not.icrilv(i) .and. taup(i,k) .gt. 0.0_r8 ) then + temv = 1.0_r8 / velco(i,k) + tem1 = coefm(i)/(dxy(i)/ncleff)*(ro(i,kp1)+ro(i,k))*brvf(i)*velco(i,k)*0.5_r8 + hd = sqrt(taup(i,k) / tem1) + fro = brvf(i) * hd * temv + ! + ! rim is the minimum-richardson number by shutts (1985) + ! + tem2 = sqrt(usqj(i,k)) + tem = 1._r8 + tem2 * fro + rim = usqj(i,k) * (1._r8-fro) / (tem * tem) + + ! + ! check stability to employ the 'saturation hypothesis' + ! of lindzen (1981) except at tropospheric downstream regions + ! + if (rim .le. ric) then ! saturation hypothesis! + if ((oa1(i) .le. 0._r8).or.(kp1 .ge. kpblmin )) then + temc = 2.0_r8 + 1.0_r8 / tem2 + hd = velco(i,k) * (2.0_r8*sqrt(temc)-temc) / brvf(i) + taup(i,kp1) = tem1 * hd * hd + ! + ! taup is restricted to monotoncally decrease + ! to avoid unexpected high taup in calculation + ! + taup(i,kp1)=min(tem1*hd*hd,taup(i,k)) + ! + ! add vertical decrease at low level below hint (Kim and Doyle 2005) + ! where Ri first decreases + ! + if (scorer_on.and.k.gt.klowtop(i).and.zl(i,k).le.zl_hint(i).and.k.lt.kte-1) then + l1=(9.81_r8*bnv2(i,kp1)/velco(i,kp1)**2) + l2=(9.81_r8*bnv2(i,k)/velco(i,k)**2) + taup(i,kp1)=min(taup(i,k),taup(i,k)*(l1/l2),tem1*hd*hd) + endif endif + else ! no wavebreaking! + taup(i,kp1) = taup(i,k) endif - else ! no wavebreaking! - taup(i,kp1) = taup(i,k) endif endif - endif + enddo enddo - enddo - if(lcap.lt.kte) then - do klcap = lcapp1,kte - do i = its,ite - taup(i,klcap) = prsi(i,klcap) / prsi(i,lcap) * taup(i,lcap) + if(lcap.lt.kte) then + do klcap = lcapp1,kte + do i = its,ite + taup(i,klcap) = prsi(i,klcap) / prsi(i,lcap) * taup(i,lcap) + enddo enddo - enddo - endif + endif ENDIF !END LARGE-SCALE TAU CALCULATION !=============================================================== @@ -1349,11 +1495,11 @@ subroutine od2d(dudt,dvdt,dthdt,ncleff,ncd,sncleff, & !=============================================================== IF (gsd_gwd_bl.and.(ls_taper .GT. 1.E-02)) THEN - do i = its,ite - if(.not.ldrag(i)) then - ! - !------- determine the height of flow-blocking layer - ! + do i = its,ite + if(.not.ldrag(i)) then + ! + !------- determine the height of flow-blocking layer + ! kblk = 0 pe = 0.0_r8 @@ -1364,9 +1510,9 @@ subroutine od2d(dudt,dvdt,dthdt,ncleff,ncd,sncleff, & !divided by g*ro is to turn del(pa) into height pe = pe + bnv2(i,k)*(zl(i,komax(i))-zl(i,k))*del(i,k)/g/ro(i,k) ke = 0.5_r8*((rcs*u1(i,k))**2._r8+(rcs*v1(i,k))**2._r8) - ! - !---------- apply flow-blocking drag when pe >= ke - ! + ! + !---------- apply flow-blocking drag when pe >= ke + ! if(pe.ge.ke) then kblk = k kblk = min(kblk,kbl(i)) @@ -1376,9 +1522,9 @@ subroutine od2d(dudt,dvdt,dthdt,ncleff,ncd,sncleff, & enddo if(kblk.ne.0) then - ! - !--------- compute flow-blocking stress - ! + ! + !--------- compute flow-blocking stress + ! !dxmax_ls is different than the usual one !because the taper is very different @@ -1403,85 +1549,82 @@ subroutine od2d(dudt,dvdt,dthdt,ncleff,ncd,sncleff, & ! !taup(i,:) = taup(i,:) + taufb(i,:) ! Keep taup and taufb separate for now endif - endif - enddo + endif + enddo ENDIF ! end blocking drag !=========================================================== IF (gsd_gwd_ls.OR.gsd_gwd_bl.and.(ls_taper .GT. 1.E-02)) THEN - ! - ! calculate - (g)*d(tau)/d(pressure) and deceleration terms dtaux, dtauy - ! - - do k = kts,kte - do i = its,ite - taud_ls(i,k) = 1._r8 * (taup(i,k+1) - taup(i,k)) * csg / del(i,k) - taud_bl(i,k) = 1._r8 * (taufb(i,k+1) - taufb(i,k)) * csg / del(i,k) - enddo - enddo - ! - ! limit de-acceleration (momentum deposition ) at top to 1/2 value - ! the idea is some stuff must go out the 'top' - ! - - do klcap = lcap,kte - do i = its,ite - taud_ls(i,klcap) = taud_ls(i,klcap) * factop - taud_bl(i,klcap) = taud_bl(i,klcap) * factop - enddo - enddo - - ! - ! if the gravity wave drag would force a critical line - ! in the lower ksmm1 layers during the next deltim timestep, - ! then only apply drag until that critical line is reached. - ! - do k = kts,kpblmax-1 + ! + ! calculate - (g)*d(tau)/d(pressure) and deceleration terms dtaux, dtauy + ! + do k = kts,kte + do i = its,ite + taud_ls(i,k) = 1._r8 * (taup(i,k+1) - taup(i,k)) * csg / del(i,k) + taud_bl(i,k) = 1._r8 * (taufb(i,k+1) - taufb(i,k)) * csg / del(i,k) + enddo + enddo + ! + ! limit de-acceleration (momentum deposition ) at top to 1/2 value + ! the idea is some stuff must go out the 'top' + ! + do klcap = lcap,kte do i = its,ite - if (k .le. kbl(i)) then - if((taud_ls(i,k)+taud_bl(i,k)).ne.0._r8) & - dtfac(i) = min(dtfac(i),abs(velco(i,k) & - /(deltim*rcs*(taud_ls(i,k)+taud_bl(i,k))))) - endif + taud_ls(i,klcap) = taud_ls(i,klcap) * factop + taud_bl(i,klcap) = taud_bl(i,klcap) * factop enddo - enddo + enddo + ! + ! if the gravity wave drag would force a critical line + ! in the lower ksmm1 layers during the next deltim timestep, + ! then only apply drag until that critical line is reached. + ! + do k = kts,kpblmax-1 + do i = its,ite + if (k .le. kbl(i)) then + if((taud_ls(i,k)+taud_bl(i,k)).ne.0._r8) & + dtfac(i) = min(dtfac(i),abs(velco(i,k) & + /(deltim*rcs*(taud_ls(i,k)+taud_bl(i,k))))) + endif + enddo + enddo - do k = kts,kte + do k = kts,kte do i = its,ite - taud_ls(i,k) = taud_ls(i,k) * dtfac(i) * ls_taper - !apply limiter for ogwd - !1.dudt < |c-u|/dt, so u-c cannot change sign(u^n+1 = u^n + du/dt * dt) - !2.dudt shr_kind_r8 - use ppgrid, only: pcols, pver, psubcols,nvar_dirOA,nvar_dirOL + use ppgrid, only: pcols, pver, psubcols use constituents, only: pcnst, qmin, cnst_name, icldliq, icldice use geopotential, only: geopotential_t use physconst, only: zvir, gravit, cpair, rair, cpairv, rairv @@ -137,16 +137,6 @@ module physics_types cid ! unique column id integer :: ulatcnt, &! number of unique lats in chunk uloncnt ! number of unique lons in chunk - real(r8), dimension(:),allocatable :: & - oc !convexity of high-res grid height - real(r8), dimension(:,:),allocatable :: & - oadir !orographic asymmetry in a coarse grid - real(r8), dimension(:,:),allocatable :: & - ol !orographic length in a coarse grid - real(r8), dimension(:),allocatable :: & - pblh !get plantet boundary layer height - real(r8), dimension(:),allocatable :: & - ribulk end type physics_state !------------------------------------------------------------------------------- @@ -1839,21 +1829,7 @@ subroutine physics_state_alloc(state,lchnk,psetcols) allocate(state%cid(psetcols), stat=ierr) if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation error for state%cid') - allocate(state%oc(psetcols), stat=ierr) - if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation error for state%oc') - allocate(state%oadir(psetcols,nvar_dirOA), stat=ierr) - if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation error for state%oadir') - allocate(state%ol(psetcols,nvar_dirOL), stat=ierr) - if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation error for state%ol') - allocate(state%pblh(psetcols), stat=ierr) - if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation error for state%pblh') - allocate(state%ribulk(psetcols), stat=ierr) - if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation error for state%ribulk') - state%oc(:)=inf - state%oadir(:,:)=inf - state%ol(:,:)=inf - state%pblh(:)=inf - state%ribulk(:)=0.0_r8!inf + state%lat(:) = inf state%lon(:) = inf state%ulat(:) = inf diff --git a/components/eam/src/physics/cam/physpkg.F90 b/components/eam/src/physics/cam/physpkg.F90 index c7b8da3c8938..1dd0e69a6d33 100644 --- a/components/eam/src/physics/cam/physpkg.F90 +++ b/components/eam/src/physics/cam/physpkg.F90 @@ -156,6 +156,7 @@ subroutine phys_register use radiation, only: radiation_register use co2_cycle, only: co2_register use co2_diagnostics, only: co2_diags_register + use gw_drag, only: gw_register use flux_avg, only: flux_avg_register use iondrag, only: iondrag_register use ionosphere, only: ionos_register @@ -316,6 +317,8 @@ subroutine phys_register call co2_register() call co2_diags_register() + call gw_register() + ! register data model ozone with pbuf if (cam3_ozone_data_on) then call cam3_ozone_data_register() @@ -906,7 +909,7 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_out ) ! CAM3 prescribed ozone if (cam3_ozone_data_on) call cam3_ozone_data_init(phys_state) - call gw_init() + call gw_init(pbuf2d) call rayleigh_friction_init() @@ -1321,7 +1324,7 @@ subroutine phys_run2(phys_state, ztodt, phys_tend, pbuf2d, cam_out, & use cam_diagnostics,only: diag_deallocate, diag_surf - use comsrf, only: trefmxav, trefmnav, sgh, sgh30, fsds, oc, oadir, ol + use comsrf, only: trefmxav, trefmnav, sgh, sgh30, fsds use physconst, only: stebol, latvap #if ( defined OFFLINE_DYN ) use metdata, only: get_met_srf2 @@ -1329,7 +1332,7 @@ subroutine phys_run2(phys_state, ztodt, phys_tend, pbuf2d, cam_out, & use time_manager, only: get_nstep, is_first_step, is_end_curr_month, & is_first_restart_step, is_last_step use check_energy, only: ieflx_gmean, check_ieflx_fix - use phys_control, only: ieflx_opt,use_od_ls,use_od_bl + use phys_control, only: ieflx_opt use co2_diagnostics,only: get_total_carbon, print_global_carbon_diags, & co2_diags_store_fields, co2_diags_read_fields use co2_cycle, only: co2_transport @@ -1432,13 +1435,7 @@ subroutine phys_run2(phys_state, ztodt, phys_tend, pbuf2d, cam_out, & call t_startf('diag_surf') call diag_surf(cam_in(c), cam_out(c), phys_state(c)%ps,trefmxav(1,c), trefmnav(1,c)) call t_stopf('diag_surf') - ! for tranport of ogwd related parameters - if ( use_od_ls .or. use_od_bl ) then - phys_state(c)%oc (:) =oc (:,c) - phys_state(c)%oadir(:,:) =oadir (:,:,c) - phys_state(c)%ol (:,:) =ol (:,:,c) - endif - ! + call tphysac(ztodt, cam_in(c), & sgh(1,c), sgh30(1,c), cam_out(c), & phys_state(c), phys_tend(c), phys_buffer_chunk, phys_diag(c), & @@ -1840,7 +1837,7 @@ subroutine tphysac (ztodt, cam_in, & ! If CLUBB is called, do not call vertical diffusion, but still ! calculate surface friction velocity (ustar) and Obukhov length - call clubb_surface ( state, cam_in, surfric, obklen) + call clubb_surface ( state, cam_in, pbuf, surfric, obklen) ! Diagnose tracer mixing ratio tendencies from surface fluxes, ! then update the mixing ratios. (If cflx_cpl_opt==2, these are done in diff --git a/components/eam/src/physics/cam/ppgrid.F90 b/components/eam/src/physics/cam/ppgrid.F90 index 8ef5d205703b..a2bbc5e7fad9 100644 --- a/components/eam/src/physics/cam/ppgrid.F90 +++ b/components/eam/src/physics/cam/ppgrid.F90 @@ -21,8 +21,6 @@ module ppgrid public psubcols public pver public pverp - public nvar_dirOA - public nvar_dirOL ! Grid point resolution parameters @@ -32,9 +30,6 @@ module ppgrid integer psubcols ! number of sub-columns (max) integer pver ! number of vertical levels integer pverp ! pver + 1 - !added for ogwd - integer nvar_dirOA - integer nvar_dirOL #ifdef PPCOLS parameter (pcols = PCOLS) @@ -42,9 +37,6 @@ module ppgrid parameter (psubcols = PSUBCOLS) parameter (pver = PLEV) parameter (pverp = pver + 1 ) - !added for ogwd - parameter (nvar_dirOA =2+1 )!avoid bug when nvar_dirOA is 2 - parameter (nvar_dirOL =180)!set for 360 degrees wind direction ! ! start, end indices for chunks owned by a given MPI task ! (set in phys_grid_init).