Add an option to calculate omega in non-hydrostatic runs similar to h…

…ydrostatic mode when 'pass_full_omega_to_physics_in_non_hydrostatic_mode' is set to true

driver/fvGFS/atmosphere.F90

@@ -192,6 +192,7 @@ module atmosphere_mod
 use fv_control_mod, only: fv_control_init, fv_end, ngrids
 use fv_eta_mod, only: get_eta_level
 use fv_fill_mod, only: fill_gfs
+use dyn_core_mod, only: del2_cubed
 use fv_dynamics_mod, only: fv_dynamics
 use fv_nesting_mod, only: twoway_nesting
 use boundary_mod, only: fill_nested_grid
@@ -709,6 +710,22 @@ subroutine atmosphere_dynamics ( Time )
  endif
 
  enddo !p_split
+ if (.not. Atm(n)%flagstruct%hydrostatic .and. .not. Atm(n)%flagstruct%pass_full_omega_to_physics_in_non_hydrostatic_mode) then
+ Atm(n)%omga(isc:iec,jsc:jec,1:npz) = Atm(n)%delp(isc:iec,jsc:jec,1:npz) / Atm(n)%delz(isc:iec,jsc:jec,1:npz) * Atm(n)%w(isc:iec,jsc:jec,1:npz)
+ if(Atm(n)%flagstruct%nf_omega>0) then
+ call del2_cubed(&
+ Atm(n)%omga, &
+ 0.18*Atm(n)%gridstruct%da_min, &
+ Atm(n)%gridstruct, &
+ Atm(n)%domain, &
+ Atm(n)%npx, &
+ Atm(n)%npy, &
+ Atm(n)%npz, &
+ Atm(n)%flagstruct%nf_omega, &
+ Atm(n)%bd)
+ endif
+ endif
+
  call mpp_clock_end (id_dynam)
 
 !-----------------------------------------------------
@@ -2052,6 +2069,7 @@ subroutine atmos_phys_driver_statein (IPD_Data, Atm_block,flip_vc)
  real(kind=kind_phys) :: pk0inv, ptop, pktop
  real(kind=kind_phys) :: rTv, dm, qgrs_rad
  integer :: nb, blen, npz, i, j, k, ix, k1, kz, dnats, nq_adv
+
 #ifdef MULTI_GASES
  real :: q_grs(nq), q_min
 #endif
@@ -2118,7 +2136,7 @@ subroutine atmos_phys_driver_statein (IPD_Data, Atm_block,flip_vc)
  if(associated(IPD_Data(nb)%Statein%wgrs) .and. .not. Atm(mygrid)%flagstruct%hydrostatic) then
  IPD_Data(nb)%Statein%wgrs(ix,k) = _DBL_(_RL_(Atm(mygrid)%w(i,j,k1)))
  endif
- IPD_Data(nb)%Statein%vvl(ix,k)  = _DBL_(_RL_(Atm(mygrid)%omga(i,j,k1)))
+ IPD_Data(nb)%Statein%vvl(ix,k) = _DBL_(_RL_(Atm(mygrid)%omga(i,j,k1)))
  IPD_Data(nb)%Statein%prsl(ix,k) = _DBL_(_RL_(Atm(mygrid)%delp(i,j,k1))) ! Total mass
  if (Atm(mygrid)%flagstruct%do_skeb)IPD_Data(nb)%Statein%diss_est(ix,k) = _DBL_(_RL_(Atm(mygrid)%diss_est(i,j,k1)))
 

model/dyn_core.F90

@@ -857,7 +857,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp,
 
  endif
 
- if( hydrostatic .and. (.not.flagstruct%use_old_omega) .and. last_step ) then
+ if( (.not.flagstruct%use_old_omega) .and. last_step ) then
 ! Average horizontal "convergence" to cell center
  do j=js,je
  do i=is,ie
@@ -894,7 +894,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp,
  nord_k, nord_v(k), nord_w, nord_t, flagstruct%dddmp, d2_divg, flagstruct%d4_bg, &
  damp_vt(k), damp_w, damp_t, d_con_k, hydrostatic, gridstruct, flagstruct, bd)
 
- if( hydrostatic .and. (.not.flagstruct%use_old_omega) .and. last_step ) then
+ if((.not.flagstruct%use_old_omega) .and. last_step ) then
 ! Average horizontal "convergence" to cell center
  do j=js,je
  do i=is,ie
@@ -1264,7 +1264,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp,
 
 #ifdef SW_DYNAMICS
 #else
- if ( hydrostatic .and. last_step ) then
+ if ( last_step ) then
  if ( flagstruct%use_old_omega ) then
 !$OMP parallel do default(none) shared(is,ie,js,je,npz,omga,pe,pem,rdt)
  do k=1,npz

model/fv_arrays.F90

@@ -926,6 +926,7 @@ module fv_arrays_mod
  integer :: nrows_blend = 0 !< # of blending rows in the outer integration domain.
  logical :: write_restart_with_bcs = .false. !< Default setting for using DA-updated BC files
  logical :: regional_bcs_from_gsi = .false. !< Default setting for writing restart files with boundary rows
+ logical :: pass_full_omega_to_physics_in_non_hydrostatic_mode = .false. !< Default to passing local omega to physics in non-hydrostatic 
 
 
  !>Convenience pointers
@@ -1510,6 +1511,7 @@ subroutine allocate_fv_atmos_type(Atm, isd_in, ied_in, jsd_in, jed_in, is_in, ie
  allocate ( Atm%ts(is:ie,js:je) )
  allocate ( Atm%phis(isd:ied ,jsd:jed ) )
  allocate ( Atm%omga(isd:ied ,jsd:jed ,npz) ); Atm%omga=0.
+
  allocate ( Atm%ua(isd:ied ,jsd:jed ,npz) )
  allocate ( Atm%va(isd:ied ,jsd:jed ,npz) )
  allocate ( Atm%uc(isd:ied+1,jsd:jed ,npz) )

model/fv_control.F90

@@ -415,6 +415,7 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split,
  logical, pointer :: write_only_coarse_intermediate_restarts
  logical, pointer :: write_coarse_agrid_vel_rst
  logical, pointer :: write_coarse_dgrid_vel_rst
+ logical, pointer :: pass_full_omega_to_physics_in_non_hydrostatic_mode
  !!!!!!!!!! END POINTERS !!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
  this_grid = -1 ! default
@@ -998,6 +999,7 @@ subroutine set_namelist_pointers(Atm)
  write_only_coarse_intermediate_restarts => Atm%coarse_graining%write_only_coarse_intermediate_restarts
  write_coarse_agrid_vel_rst => Atm%coarse_graining%write_coarse_agrid_vel_rst
  write_coarse_dgrid_vel_rst => Atm%coarse_graining%write_coarse_dgrid_vel_rst
+ pass_full_omega_to_physics_in_non_hydrostatic_mode => Atm%flagstruct%pass_full_omega_to_physics_in_non_hydrostatic_mode
  end subroutine set_namelist_pointers
 
 
@@ -1091,7 +1093,8 @@ subroutine read_namelist_fv_core_nml(Atm)
  write_coarse_restart_files,&
  write_coarse_diagnostics,&
  write_only_coarse_intermediate_restarts, &
- write_coarse_agrid_vel_rst, write_coarse_dgrid_vel_rst, ignore_rst_cksum
+ write_coarse_agrid_vel_rst, write_coarse_dgrid_vel_rst, &
+ pass_full_omega_to_physics_in_non_hydrostatic_mode, ignore_rst_cksum
 
 
  ! Read FVCORE namelist

model/fv_dynamics.F90

@@ -293,7 +293,7 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill,
  integer :: sphum, liq_wat = -999, ice_wat = -999 ! GFDL physics
  integer :: rainwat = -999, snowwat = -999, graupel = -999, hailwat = -999, cld_amt = -999
  integer :: theta_d = -999
- logical used, do_omega
+ logical used
  integer, parameter :: max_packs=13
  type(group_halo_update_type), save :: i_pack(max_packs)
  integer :: is, ie, js, je
@@ -763,7 +763,6 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill,
  if ( iq==cld_amt ) kord_tracer(iq) = 9 ! monotonic
  enddo
 
- do_omega = hydrostatic .and. last_step
  call timing_on('Remapping')
 #ifdef AVEC_TIMERS
  call avec_timer_start(6)
@@ -777,7 +776,7 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill,
  ng, ua, va, omga, dp1, ws, fill, reproduce_sum, &
  idiag%id_mdt>0, dtdt_m, ptop, ak, bk, pfull, gridstruct, domain, &
  flagstruct%do_sat_adj, hydrostatic, flagstruct%phys_hydrostatic, &
- hybrid_z, do_omega,  &
+ hybrid_z, &
  flagstruct%adiabatic, do_adiabatic_init, flagstruct%do_inline_mp, &
  inline_mp, flagstruct%c2l_ord, bd, flagstruct%fv_debug, &
  flagstruct%moist_phys)
@@ -820,22 +819,13 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill,
  endif
 #endif
 
- if( last_step ) then
- if( .not. hydrostatic ) then
-!$OMP parallel do default(none) shared(is,ie,js,je,npz,omga,delp,delz,w)
- do k=1,npz
- do j=js,je
- do i=is,ie
- omga(i,j,k) = delp(i,j,k)/delz(i,j,k)*w(i,j,k)
- enddo
- enddo
- enddo
- endif
 !--------------------------
 ! Filter omega for physics:
 !--------------------------
- if(flagstruct%nf_omega>0) &
- call del2_cubed(omga, 0.18*gridstruct%da_min, gridstruct, domain, npx, npy, npz, flagstruct%nf_omega, bd)
+ if( last_step ) then
+ if(flagstruct%nf_omega>0) then 
+ call del2_cubed(omga, 0.18*gridstruct%da_min, gridstruct, domain, npx, npy, npz, flagstruct%nf_omega, bd)
+ endif
  endif
  end if
 #endif

model/fv_mapz.F90

@@ -137,7 +137,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, &
  akap, cappa, kord_mt, kord_wz, kord_tr, kord_tm, peln, te0_2d, &
  ng, ua, va, omga, te, ws, fill, reproduce_sum, out_dt, dtdt, &
  ptop, ak, bk, pfull, gridstruct, domain, do_sat_adj, &
- hydrostatic, phys_hydrostatic, hybrid_z, do_omega, adiabatic, do_adiabatic_init, &
+ hydrostatic, phys_hydrostatic, hybrid_z, adiabatic, do_adiabatic_init, &
  do_inline_mp, inline_mp, c2l_ord, bd, fv_debug, &
  moist_phys)
  logical, intent(in):: last_step
@@ -169,7 +169,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, &
  logical, intent(in):: do_inline_mp
  logical, intent(in):: fill !< fill negative tracers
  logical, intent(in):: reproduce_sum
- logical, intent(in):: do_omega, adiabatic, do_adiabatic_init
+ logical, intent(in):: adiabatic, do_adiabatic_init
  real, intent(in) :: ptop
  real, intent(in) :: ak(km+1)
  real, intent(in) :: bk(km+1)
@@ -258,7 +258,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, &
 #ifdef MULTI_GASES
 !$OMP num_gas, &
 #endif
-!$OMP hs,w,ws,kord_wz,do_omega,omga,rrg,kord_mt,pe4) &
+!$OMP hs,w,ws,kord_wz,omga,rrg,kord_mt,pe4) &
 !$OMP private(qv,gz,cvm,kp,k_next,bkh,dp2, &
 !$OMP pe0,pe1,pe2,pe3,pk1,pk2,pn2,phis,q2,w2)
  do 1000 j=js,je+1
@@ -508,8 +508,8 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, &
  enddo
 
 !----------------
- if ( do_omega ) then
-! Start do_omega
+ if ( last_step ) then
+! Start last_step
 ! Copy omega field to pe3
  do i=is,ie
  pe3(i,1) = 0.
@@ -588,7 +588,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, &
  endif
 
 ! Interpolate omega/pe3 (defined at pe0) to remapped cell center (dp2)
- if ( do_omega ) then
+ if ( last_step ) then
  do k=1,km
  do i=is,ie
  dp2(i,k) = 0.5*(peln(i,k,j) + peln(i,k+1,j))
@@ -608,7 +608,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, &
  enddo
  enddo
  enddo
- endif ! end do_omega
+ endif ! end last_step
 
  endif !(j < je+1)