diff --git a/driver/UFS/atmosphere.F90 b/driver/UFS/atmosphere.F90 index b9abdb28d..61eb52d28 100644 --- a/driver/UFS/atmosphere.F90 +++ b/driver/UFS/atmosphere.F90 @@ -182,6 +182,8 @@ module atmosphere_mod #ifdef GFS_TYPES use GFS_typedefs, only: IPD_control_type => GFS_control_type, kind_phys use GFS_typedefs, only: GFS_statein_type, GFS_stateout_type, GFS_sfcprop_type +! SA-3D-TKE (kyf) +use GFS_typedefs, only: GFS_tbd_type #else use IPD_typedefs, only: IPD_data_type, IPD_control_type, kind_phys => IPD_kind_phys #endif @@ -697,7 +699,10 @@ subroutine atmosphere_dynamics ( Time ) Atm(n)%pe, Atm(n)%pk, Atm(n)%peln, & Atm(n)%pkz, Atm(n)%phis, Atm(n)%q_con, & Atm(n)%omga, Atm(n)%ua, Atm(n)%va, Atm(n)%uc, & - Atm(n)%vc, Atm(n)%ak, Atm(n)%bk, Atm(n)%mfx, & + Atm(n)%vc, & +!The following variable is used for SA-3D-TKE (kyf) (modify for data structure) + Atm(n)%sa3dtke_var, & + Atm(n)%ak, Atm(n)%bk, Atm(n)%mfx, & Atm(n)%mfy , Atm(n)%cx, Atm(n)%cy, Atm(n)%ze0, & Atm(n)%flagstruct%hybrid_z, & Atm(n)%gridstruct, Atm(n)%flagstruct, & @@ -1897,6 +1902,8 @@ subroutine adiabatic_init(zvir,nudge_dz,time) Atm(mygrid)%pt, Atm(mygrid)%delp, Atm(mygrid)%q, Atm(mygrid)%ps, & Atm(mygrid)%pe, Atm(mygrid)%pk, Atm(mygrid)%peln, Atm(mygrid)%pkz, Atm(mygrid)%phis, & Atm(mygrid)%q_con, Atm(mygrid)%omga, Atm(mygrid)%ua, Atm(mygrid)%va, Atm(mygrid)%uc, Atm(mygrid)%vc, & +!The following variable is used for SA-3D-TKE (kyf) (modify for data structure) + Atm(mygrid)%sa3dtke_var, & Atm(mygrid)%ak, Atm(mygrid)%bk, Atm(mygrid)%mfx, Atm(mygrid)%mfy, & Atm(mygrid)%cx, Atm(mygrid)%cy, Atm(mygrid)%ze0, Atm(mygrid)%flagstruct%hybrid_z, & Atm(mygrid)%gridstruct, Atm(mygrid)%flagstruct, & @@ -1912,6 +1919,8 @@ subroutine adiabatic_init(zvir,nudge_dz,time) Atm(mygrid)%pt, Atm(mygrid)%delp, Atm(mygrid)%q, Atm(mygrid)%ps, & Atm(mygrid)%pe, Atm(mygrid)%pk, Atm(mygrid)%peln, Atm(mygrid)%pkz, Atm(mygrid)%phis, & Atm(mygrid)%q_con, Atm(mygrid)%omga, Atm(mygrid)%ua, Atm(mygrid)%va, Atm(mygrid)%uc, Atm(mygrid)%vc, & +!The following three variables are used for SA-3D-TKE (kyf) (modify for data structure) + Atm(mygrid)%sa3dtke_var, & Atm(mygrid)%ak, Atm(mygrid)%bk, Atm(mygrid)%mfx, Atm(mygrid)%mfy, & Atm(mygrid)%cx, Atm(mygrid)%cy, Atm(mygrid)%ze0, Atm(mygrid)%flagstruct%hybrid_z, & Atm(mygrid)%gridstruct, Atm(mygrid)%flagstruct, & @@ -1988,6 +1997,8 @@ subroutine adiabatic_init(zvir,nudge_dz,time) Atm(mygrid)%pt, Atm(mygrid)%delp, Atm(mygrid)%q, Atm(mygrid)%ps, & Atm(mygrid)%pe, Atm(mygrid)%pk, Atm(mygrid)%peln, Atm(mygrid)%pkz, Atm(mygrid)%phis, & Atm(mygrid)%q_con, Atm(mygrid)%omga, Atm(mygrid)%ua, Atm(mygrid)%va, Atm(mygrid)%uc, Atm(mygrid)%vc, & +!The following three variables are used for SA-3D-TKE (kyf) (modify for data structure) + Atm(mygrid)%sa3dtke_var, & Atm(mygrid)%ak, Atm(mygrid)%bk, Atm(mygrid)%mfx, Atm(mygrid)%mfy, & Atm(mygrid)%cx, Atm(mygrid)%cy, Atm(mygrid)%ze0, Atm(mygrid)%flagstruct%hybrid_z, & Atm(mygrid)%gridstruct, Atm(mygrid)%flagstruct, & @@ -2002,6 +2013,8 @@ subroutine adiabatic_init(zvir,nudge_dz,time) Atm(mygrid)%pt, Atm(mygrid)%delp, Atm(mygrid)%q, Atm(mygrid)%ps, & Atm(mygrid)%pe, Atm(mygrid)%pk, Atm(mygrid)%peln, Atm(mygrid)%pkz, Atm(mygrid)%phis, & Atm(mygrid)%q_con, Atm(mygrid)%omga, Atm(mygrid)%ua, Atm(mygrid)%va, Atm(mygrid)%uc, Atm(mygrid)%vc, & +!The following three variables are used for SA-3D-TKE (kyf) (modify for data structure) + Atm(mygrid)%sa3dtke_var, & Atm(mygrid)%ak, Atm(mygrid)%bk, Atm(mygrid)%mfx, Atm(mygrid)%mfy, & Atm(mygrid)%cx, Atm(mygrid)%cy, Atm(mygrid)%ze0, Atm(mygrid)%flagstruct%hybrid_z, & Atm(mygrid)%gridstruct, Atm(mygrid)%flagstruct, & @@ -2064,6 +2077,7 @@ end subroutine adiabatic_init !>@detail Performs a mass adjustment to be consistent with the !! GFS physics and if necessary, converts quantities to hydrostatic !! representation. +!! SA-3D-TKE (added IPD_Tbd in the input) (kyf) #if defined(OVERLOAD_R4) #define _DBL_(X) DBLE(X) #define _RL_(X) REAL(X,KIND=4) @@ -2071,9 +2085,11 @@ end subroutine adiabatic_init #define _DBL_(X) X #define _RL_(X) X #endif - subroutine atmos_phys_driver_statein (IPD_Control, IPD_Statein, Atm_block,flip_vc) + subroutine atmos_phys_driver_statein (IPD_Control, IPD_Statein, IPD_Tbd, Atm_block,flip_vc) type (IPD_control_type), intent(in) :: IPD_Control type (GFS_statein_type), intent(inout) :: IPD_Statein + ! SA-3D-TKE (added IPD_Tbd) (kyf) + type (GFS_tbd_type), intent(inout) :: IPD_Tbd type (block_control_type), intent(in) :: Atm_block logical, intent(in) :: flip_vc !-------------------------------------- @@ -2113,8 +2129,10 @@ subroutine atmos_phys_driver_statein (IPD_Control, IPD_Statein, Atm_block,flip_v !--------------------------------------------------------------------- ! use most up to date atmospheric properties when running serially !--------------------------------------------------------------------- +! SA-3D-TKE added IPD_Tbd (kyf) !$OMP parallel do default (none) & !$OMP shared (Atm_block, Atm, IPD_Control, IPD_Statein, npz, nq, ncnst, sphum, liq_wat, & +!$OMP IPD_Tbd, & !$OMP ice_wat, rainwat, snowwat, graupel, pk0inv, ptop, & !$OMP pktop, zvir, mygrid, dnats, nq_adv, flip_vc) & #ifdef MULTI_GASES @@ -2130,6 +2148,27 @@ subroutine atmos_phys_driver_statein (IPD_Control, IPD_Statein, Atm_block,flip_v ! log(pe) <-- prsik blen = Atm_block%blksz(nb) +!The following is for SA-3D-TKE + if(Atm(mygrid)%flagstruct%sa3dtke_dyco) then +!The following is for SA-3D-TKE (pass dku to dyn_core) + do k = 1, npz + kz = npz+1-k + if(flip_vc) then + k1 = kz ! flipping the index + else + k1 = k + endif + do ix = 1, blen + i = Atm_block%index(nb)%ii(ix) + j = Atm_block%index(nb)%jj(ix) + ! SA-3D-TKE (added im) (kyf) + im = IPD_control%chunk_begin(nb)+ix-1 + ! SA-3D-TKE (modified IPD_Tbd and ix into im) (kyf) (modify for data structure) + Atm(mygrid)%sa3dtke_var%dku3d_h(i,j,k) = IPD_Tbd%dku3d_h(im,k1) + Atm(mygrid)%sa3dtke_var%dku3d_e(i,j,k) = IPD_Tbd%dku3d_e(im,k1) + enddo + enddo + endif do k = 1, npz !Indices for FV's vertical coordinate, for which 1 = top @@ -2154,6 +2193,13 @@ subroutine atmos_phys_driver_statein (IPD_Control, IPD_Statein, Atm_block,flip_v endif IPD_Statein%vvl(im,k) = _DBL_(_RL_(Atm(mygrid)%omga(i,j,k1))) IPD_Statein%prsl(im,k) = _DBL_(_RL_(Atm(mygrid)%delp(i,j,k1))) ! Total mass +!The following is for SA-3D-TKE (kyf) (modify for data structure) + if(Atm(mygrid)%flagstruct%sa3dtke_dyco) then + IPD_Statein%def_1(ix,k) = _DBL_(_RL_(Atm(mygrid)%sa3dtke_var%deform_1(i,j,k1))) + IPD_Statein%def_2(ix,k) = _DBL_(_RL_(Atm(mygrid)%sa3dtke_var%deform_2(i,j,k1))) + IPD_Statein%def_3(ix,k) = _DBL_(_RL_(Atm(mygrid)%sa3dtke_var%deform_3(i,j,k1))) + endif + if (Atm(mygrid)%flagstruct%do_skeb) IPD_Statein%diss_est(im,k) = _DBL_(_RL_(Atm(mygrid)%diss_est(i,j,k1))) if(flip_vc) then diff --git a/model/dyn_core.F90 b/model/dyn_core.F90 index 7f69a776c..d1d60d106 100644 --- a/model/dyn_core.F90 +++ b/model/dyn_core.F90 @@ -116,6 +116,7 @@ module dyn_core_mod use sw_core_mod, only: c_sw, d_sw, d_md use a2b_edge_mod, only: a2b_ord2, a2b_ord4 use nh_core_mod, only: Riem_Solver3, Riem_Solver_C, update_dz_c, update_dz_d, nh_bc + use nh_utils_mod, only: edge_profile1 ! KGao: for dudz,dvdz,dwdz calculations use tp_core_mod, only: copy_corners use fv_timing_mod, only: timing_on, timing_off use fv_diagnostics_mod, only: prt_maxmin, fv_time, prt_mxm @@ -131,11 +132,15 @@ module dyn_core_mod use diag_manager_mod, only: send_data use fv_arrays_mod, only: fv_grid_type, fv_flags_type, fv_nest_type, fv_diag_type, & fv_grid_bounds_type, R_GRID, fv_nest_BC_type_3d + use fv_arrays_mod, only: sa3dtke_type ! for SA-3D-TKE (kyf) (modify for data structure) use boundary_mod, only: extrapolation_BC, nested_grid_BC_apply_intT use fv_regional_mod, only: regional_boundary_update use fv_regional_mod, only: current_time_in_seconds, bc_time_interval use fv_regional_mod, only: delz_regBC ! TEMPORARY --- lmh +!The following is for SA-3D-TKE + use tracer_manager_mod, only: get_tracer_names, get_number_tracers, get_tracer_index + use field_manager_mod, only: MODEL_ATMOS #ifdef SW_DYNAMICS use test_cases_mod, only: test_case, case9_forcing1, case9_forcing2 @@ -178,7 +183,10 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, #endif grav, hydrostatic, & u, v, w, delz, pt, q, delp, pe, pk, phis, ws, omga, ptop, pfull, ua, va, & - uc, vc, mfx, mfy, cx, cy, pkz, peln, q_con, ak, bk, & + uc, vc, & +!The following variable is for SA-3D-TKE (kyf) (modify for data structure) + sa3dtke_var, & + mfx, mfy, cx, cy, pkz, peln, q_con, ak, bk, & ks, gridstruct, flagstruct, neststruct, idiag, bd, domain, & init_step, i_pack, end_step, diss_est,time_total) @@ -235,6 +243,9 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, real, intent(inout):: uc(bd%isd:bd%ied+1,bd%jsd:bd%jed ,npz) !< (uc, vc) are mostly used as the C grid winds real, intent(inout):: vc(bd%isd:bd%ied ,bd%jsd:bd%jed+1,npz) real, intent(inout), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: ua, va +!The following 2 variables are for SA-3D-TKE (kyf) (modify for data structure) + type(sa3dtke_type), intent(inout) :: sa3dtke_var + real, intent(inout):: q_con(bd%isd:, bd%jsd:, 1:) ! The Flux capacitors: accumulated Mass flux arrays @@ -288,6 +299,9 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, logical used real :: split_timestep_bc +!The following is for SA-3D-TKE + integer :: sgs_tke + integer :: is, ie, js, je integer :: isd, ied, jsd, jed @@ -765,13 +779,17 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, endif +!The following is modified for SA-3D-TKE +!-- add sa3dtke_var in parallel calculation +! replce dku3d_h by sa3dtke_var (kyf) (modify for data structure) call timing_on('d_sw') !$OMP parallel do default(none) shared(npz,flagstruct,nord_v,pfull,damp_vt,hydrostatic,last_step, & !$OMP is,ie,js,je,isd,ied,jsd,jed,omga,delp,gridstruct,npx,npy, & !$OMP ng,zh,vt,ptc,pt,u,v,w,uc,vc,ua,va,divgd,mfx,mfy,cx,cy, & +!$OMP sa3dtke_var, & !$OMP crx,cry,xfx,yfx,q_con,zvir,sphum,nq,q,dt,bd,rdt,iep1,jep1, & -!$OMP heat_source,diss_est,ptop,first_call) & +!$OMP heat_source,diss_est,ptop,first_call) & !$OMP private(nord_k, nord_w, nord_t, damp_w, damp_t, d2_divg, & !$OMP d_con_k,kgb, hord_m, hord_v, hord_t, hord_p, wk, heat_s,diss_e, z_rat) do k=1,npz @@ -879,20 +897,40 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, enddo enddo endif - call d_sw(vt(isd,jsd,k), delp(isd,jsd,k), ptc(isd,jsd,k), pt(isd,jsd,k), & - u(isd,jsd,k), v(isd,jsd,k), w(isd:,jsd:,k), uc(isd,jsd,k), & - vc(isd,jsd,k), ua(isd,jsd,k), va(isd,jsd,k), divgd(isd,jsd,k), & - mfx(is, js, k), mfy(is, js, k), cx(is, jsd,k), cy(isd,js, k), & - crx(is, jsd,k), cry(isd,js, k), xfx(is, jsd,k), yfx(isd,js, k), & + + if(flagstruct%sa3dtke_dyco) then + call d_sw(vt(isd,jsd,k), delp(isd,jsd,k), ptc(isd,jsd,k), pt(isd,jsd,k), & + u(isd,jsd,k), v(isd,jsd,k), w(isd:,jsd:,k), uc(isd,jsd,k), & + vc(isd,jsd,k), ua(isd,jsd,k), va(isd,jsd,k), divgd(isd,jsd,k), & + mfx(is, js, k), mfy(is, js, k), cx(is, jsd,k), cy(isd,js, k), & + crx(is, jsd,k), cry(isd,js, k), xfx(is, jsd,k), yfx(isd,js, k), & +#ifdef USE_COND + q_con(isd:,jsd:,k), z_rat(isd,jsd), & +#else + q_con(isd:,jsd:,1), z_rat(isd,jsd), & +#endif + kgb, heat_s, diss_e,zvir, sphum, nq, q, k, npz, flagstruct%inline_q, dt, & + flagstruct%hord_tr, hord_m, hord_v, hord_t, hord_p, & + 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, & +!The following optional variable is for SA-3D-TKE (kyf) + sa3dtke_var%dku3d_h(isd, jsd, k) ) + else + call d_sw(vt(isd,jsd,k), delp(isd,jsd,k), ptc(isd,jsd,k), pt(isd,jsd,k), & + u(isd,jsd,k), v(isd,jsd,k), w(isd:,jsd:,k), uc(isd,jsd,k), & + vc(isd,jsd,k), ua(isd,jsd,k), va(isd,jsd,k), divgd(isd,jsd,k), & + mfx(is, js, k), mfy(is, js, k), cx(is, jsd,k), cy(isd,js, k), & + crx(is, jsd,k), cry(isd,js, k), xfx(is, jsd,k), yfx(isd,js, k), & #ifdef USE_COND - q_con(isd:,jsd:,k), z_rat(isd,jsd), & + q_con(isd:,jsd:,k), z_rat(isd,jsd), & #else - q_con(isd:,jsd:,1), z_rat(isd,jsd), & + q_con(isd:,jsd:,1), z_rat(isd,jsd), & #endif - kgb, heat_s, diss_e,zvir, sphum, nq, q, k, npz, flagstruct%inline_q, dt, & - flagstruct%hord_tr, hord_m, hord_v, hord_t, hord_p, & - 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) + kgb, heat_s, diss_e,zvir, sphum, nq, q, k, npz, flagstruct%inline_q, dt, & + flagstruct%hord_tr, hord_m, hord_v, hord_t, hord_p, & + 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) + endif if((.not.flagstruct%use_old_omega) .and. last_step ) then ! Average horizontal "convergence" to cell center @@ -1370,6 +1408,25 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, enddo ! time split loop first_call=.false. !----------------------------------------------------- +!The following is for SA-3D-TKE (kyf) (modify for data structure)(update mpi use) +!--calculating shear deformation and TKE transport for 3d TKE scheme + if(flagstruct%sa3dtke_dyco) then + if ( end_step ) then + sgs_tke = get_tracer_index(MODEL_ATMOS, 'sgs_tke') + call mpp_update_domains(q(:,:,:,sgs_tke), domain, complete=.false.) + call mpp_update_domains(sa3dtke_var%dku3d_e(:,:,:), domain, complete=.true.) ! update dku3d_e at halo + + call diff3d(npx, npy, npz, nq, ua, va, w, & + q, sa3dtke_var%deform_1, sa3dtke_var%deform_2, & + sa3dtke_var%deform_3, sa3dtke_var%dku3d_e, & + zh, dp_ref, gridstruct, bd) + call mpp_update_domains(sa3dtke_var%deform_1, domain, complete=.false.) + call mpp_update_domains(sa3dtke_var%deform_2, domain, complete=.false.) + call mpp_update_domains(sa3dtke_var%deform_3, domain, complete=.true.) + endif + endif + + if ( nq > 0 .and. .not. flagstruct%inline_q ) then call timing_on('COMM_TOTAL') call timing_on('COMM_TRACER') @@ -2112,6 +2169,286 @@ subroutine one_grad_p(u, v, pk, gz, divg2, delp, dt, ng, gridstruct, bd, npx, np end subroutine one_grad_p +!The following is for SA-3D-TKE +subroutine diff3d(npx, npy, npz, nq, ua, va, w, q, & + deform_1, deform_2, deform_3, dku3d_e, & + zh, dp_ref, gridstruct, bd) + + use fv_arrays_mod, only: fv_grid_type + + integer, intent(in) :: nq, npx, npy, npz + type(fv_grid_bounds_type), intent(IN) :: bd + real, intent(in) :: ua(bd%isd:bd%ied, bd%jsd:bd%jed, npz) + real, intent(in) :: va(bd%isd:bd%ied, bd%jsd:bd%jed, npz) + real, intent(in) :: w(bd%isd:bd%ied, bd%jsd:bd%jed, npz) + real, intent(in) :: q(bd%isd:bd%ied, bd%jsd:bd%jed, npz, nq) + real, intent(in) :: zh(bd%isd:bd%ied, bd%jsd:bd%jed, npz+1) + real, intent(in) :: dku3d_e(bd%isd:bd%ied ,bd%jsd:bd%jed,npz) + + !real, intent(in) :: gz(bd%is:,bd%js:,1:) ! KGao: dims may not be right; zh is now used + real, intent(in) :: dp_ref(npz) + + real, intent(out) :: deform_1(bd%isd:bd%ied, bd%jsd:bd%jed, npz) + real, intent(out) :: deform_2(bd%isd:bd%ied, bd%jsd:bd%jed, npz) + real, intent(out) :: deform_3(bd%isd:bd%ied, bd%jsd:bd%jed, npz) + type(fv_grid_type), intent(IN), target :: gridstruct + +! local + + real, pointer, dimension(:,:) :: dx, dy, rarea + integer :: is, ie, js, je, isd, ied, jsd, jed, i, j, k + real:: dudx(bd%isd:bd%ied,bd%jsd:bd%jed,npz) + real:: dudy(bd%isd:bd%ied,bd%jsd:bd%jed,npz) + real:: dvdx(bd%isd:bd%ied,bd%jsd:bd%jed,npz) + real:: dvdy(bd%isd:bd%ied,bd%jsd:bd%jed,npz) + real:: dwdx(bd%isd:bd%ied,bd%jsd:bd%jed,npz) + real:: dwdy(bd%isd:bd%ied,bd%jsd:bd%jed,npz) + real:: dudz(bd%isd:bd%ied,bd%jsd:bd%jed,npz) + real:: dvdz(bd%isd:bd%ied,bd%jsd:bd%jed,npz) + real:: dwdz(bd%isd:bd%ied,bd%jsd:bd%jed,npz) + real:: ut(bd%isd:bd%ied+1,bd%jsd:bd%jed) + real:: vt(bd%isd:bd%ied, bd%jsd:bd%jed+1) + real:: u_e(bd%is:bd%ie,npz+1) + real:: v_e(bd%is:bd%ie,npz+1) + real:: w_e(bd%is:bd%ie,npz+1) + + real:: dkdx(bd%isd:bd%ied,bd%jsd:bd%jed,npz) + real:: dkdy(bd%isd:bd%ied,bd%jsd:bd%jed,npz) + + integer :: sgs_tke + real :: tke_1(bd%isd:bd%ied+2, bd%jsd:bd%jed+2) + real :: dedy_1(bd%isd:bd%ied+1, bd%jsd:bd%jed+1, npz) + real :: dedy_2(bd%isd:bd%ied, bd%jsd:bd%jed, npz) + real :: dedx_1(bd%isd:bd%ied+1, bd%jsd:bd%jed+1, npz) + real :: dedx_2(bd%isd:bd%ied, bd%jsd:bd%jed, npz) + real :: tke_2(npz) + real :: dedz_1(npz) + real :: dedz_2(bd%isd:bd%ied, bd%jsd:bd%jed,npz) + real :: tmp + + real :: dz + + dx => gridstruct%dx + dy => gridstruct%dy + rarea => gridstruct%rarea + + is = bd%is + ie = bd%ie + js = bd%js + je = bd%je + isd = bd%isd + ied = bd%ied + jsd = bd%jsd + jed = bd%jed + +!=========================================================== +! Calculate deform_1 and deform_2 +! deform_1 = 2*dw/dz**2+(du/dz**2+dv/dz**2) +! +(dw/dx*du/dz+dw/dy*dv/dz) +! deform_2 = 2*(du/dx**2+dv/dy**2)+(du/dy+dv/dx)**2 +! +(dw/dx**2+dw/dy**2)+(dw/dx*du/dz+dw/dy*dv/dz) +!=========================================================== + +! KGao: make ut and vt private as suggested by Lucas + +!$OMP parallel do default(none) shared(npz,is,ie,js,je,ua,va,w,dx,dy,rarea, & +!$OMP dudx,dudy,dvdx,dvdy,dwdx,dwdy) & +!$OMP private(ut,vt) + do k=1,npz +!------------------------------------- +! get du/dy and dv/dx + do j=js,je+1 + do i=is,ie + vt(i,j) = ua(i,j,k)*dx(i,j) + enddo + enddo + do j=js,je + do i=is,ie+1 + ut(i,j) = va(i,j,k)*dy(i,j) + enddo + enddo + do j=js,je + do i=is,ie + dudy(i,j,k) = rarea(i,j)*(vt(i,j+1)-vt(i,j)) + dvdx(i,j,k) = rarea(i,j)*(ut(i+1,j)-ut(i,j)) + enddo + enddo +!------------------------------------- +! get du/dx and dv/dy + do j=js,je + do i=is,ie+1 + ut(i,j) = ua(i,j,k)*dy(i,j) + enddo + enddo + do j=js,je+1 + do i=is,ie + vt(i,j) = va(i,j,k)*dx(i,j) + enddo + enddo + do j=js,je + do i=is,ie + dudx(i,j,k) = rarea(i,j)*(ut(i+1,j)-ut(i,j)) + dvdy(i,j,k) = rarea(i,j)*(vt(i,j+1)-vt(i,j)) + enddo + enddo +!------------------------------------- +! get dw/dx and dw/dy + do j=js,je + do i=is,ie+1 + ut(i,j) = w(i,j,k)*dy(i,j) + enddo + enddo + do j=js,je+1 + do i=is,ie + vt(i,j) = w(i,j,k)*dx(i,j) + enddo + enddo + do j=js,je + do i=is,ie + dwdx(i,j,k) = rarea(i,j)*(ut(i+1,j)-ut(i,j)) + dwdy(i,j,k) = rarea(i,j)*(vt(i,j+1)-vt(i,j)) + enddo + enddo + enddo !z loop + +!------------------------------------- +! get du/dz, dv/dz and dw/dz + +! KGao: use Lucas's method as in compute_dudz() +!$OMP parallel do default(none) shared(npz,is,ie,js,je,ua,va,w, & +!$OMP zh,dp_ref,dudz,dvdz,dwdz) & +!$OMP private(u_e,v_e,w_e,dz) + do j=js,je + call edge_profile1(ua(is:ie,j,:), u_e, is, ie, npz, dp_ref, 1) + call edge_profile1(va(is:ie,j,:), v_e, is, ie, npz, dp_ref, 1) + call edge_profile1(w(is:ie,j,:), w_e, is, ie, npz, dp_ref, 1) + do k=1,npz + do i=is,ie + dz = zh(i,j,k) - zh(i,j,k+1) + dudz(i,j,k) = (u_e(i,k)-u_e(i,k+1))/dz + dvdz(i,j,k) = (v_e(i,k)-v_e(i,k+1))/dz + dwdz(i,j,k) = (w_e(i,k)-w_e(i,k+1))/dz + enddo + enddo + enddo + +!------------------------------------- +! get deform_1 and deform_2 based on all terms + +!$OMP parallel do default(none) shared(npz,is,ie,js,je,deform_1,deform_2, & +!$OMP dudx,dudy,dudz,dvdx,dvdy,dvdz,dwdx,dwdy,dwdz) & +!$OMP private(tmp) + do k=1,npz + do j=js,je + do i=is,ie + tmp=dwdx(i,j,k)*dudz(i,j,k)+dwdy(i,j,k)*dvdz(i,j,k) + deform_1(i,j,k)= 2*dwdz(i,j,k)**2+ & + dudz(i,j,k)**2+dvdz(i,j,k)**2+tmp + deform_2(i,j,k)= & + 2*(dudx(i,j,k)**2+dvdy(i,j,k)**2)+ & + (dudy(i,j,k)+dvdx(i,j,k))**2+ & + dwdx(i,j,k)**2+dwdy(i,j,k)**2+tmp + enddo + enddo + enddo + +!=========================================================== +! Calculate deform_3 +! deform_3 = d(kqde/dx)/dx + d(kqde/dy)/dy +! where e is tke +!=========================================================== + + sgs_tke = get_tracer_index(MODEL_ATMOS, 'sgs_tke') + +! KGao: make the 2d temporay arrays private - as suggested by Lucas + +!$OMP parallel do default(none) shared(npz,is,ie,js,je,ua,va,w,q,dx,dy,rarea, & +!$OMP dku3d_e,dkdx,dkdy,dedx_1,dedy_1,sgs_tke,dedy_2,dedx_2) & +!$OMP private(ut,vt,tke_1) + do k=1,npz +!------------------------------------- + do j=js,je+2 + do i=is,ie+2 + tke_1(i,j)=q(i,j,k,sgs_tke) + enddo + enddo + do j=js,je+2 + do i=is,ie + vt(i,j)=tke_1(i,j)*dx(i,j) + enddo + enddo + do j=js,je+1 + do i=is,ie + dedy_1(i,j,k)=rarea(i,j)*(vt(i,j+1)-vt(i,j)) + enddo + enddo + do j=js,je+1 + do i=is,ie + vt(i,j)=dedy_1(i,j,k)*dx(i,j) + enddo + enddo + do j=js,je + do i=is,ie + dedy_2(i,j,k)=rarea(i,j)*(vt(i,j+1)-vt(i,j)) + enddo + enddo +!------------------------------------- + do j=js,je + do i=is,ie+2 + ut(i,j)=tke_1(i,j)*dy(i,j) + enddo + enddo + do j=js,je + do i=is,ie+1 + dedx_1(i,j,k)=rarea(i,j)*(ut(i+1,j)-ut(i,j)) + enddo + enddo + do j=js,je + do i=is,ie+1 + ut(i,j)=dedx_1(i,j,k)*dy(i,j) + enddo + enddo + do j=js,je + do i=is,ie + dedx_2(i,j,k)=rarea(i,j)*(ut(i+1,j)-ut(i,j)) + enddo + enddo + do j=js,je+1 + do i=is,ie + vt(i,j)=dku3d_e(i,j,k)*dx(i,j) + enddo + enddo + do j=js,je + do i=is,ie + dkdy(i,j,k)=rarea(i,j)*(vt(i,j+1)-vt(i,j)) + enddo + enddo + do j=js,je + do i=is,ie+1 + ut(i,j)=dku3d_e(i,j,k)*dy(i,j) + enddo + enddo + do j=js,je + do i=is,ie + dkdx(i,j,k)=rarea(i,j)*(ut(i+1,j)-ut(i,j)) + enddo + enddo + enddo ! z loop +!------------------------------------- +! get deform_3 based on all terms + +!$OMP parallel do default(none) shared(npz,is,ie,js,je,deform_3, & +!$OMP dku3d_e,dkdx,dkdy,dedx_1,dedy_1,sgs_tke,dedx_2,dedy_2) + do k=1,npz + do j=js,je + do i=is,ie + deform_3(i,j,k)=dku3d_e(i,j,k)*(dedx_2(i,j,k)+dedy_2(i,j,k)) & + +dkdx(i,j,k)*dedx_1(i,j,k)+dkdy(i,j,k)*dedy_1(i,j,k) + enddo + enddo + enddo + +end subroutine diff3d subroutine grad1_p_update(divg2, u, v, pk, gz, dt, ng, gridstruct, bd, npx, npy, npz, ptop, beta, a2b_ord) diff --git a/model/fv_arrays.F90 b/model/fv_arrays.F90 index 03a04035c..15f4daa8c 100644 --- a/model/fv_arrays.F90 +++ b/model/fv_arrays.F90 @@ -931,6 +931,8 @@ module fv_arrays_mod 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 +!This logical variable is used for SA-3D-TKE + logical :: sa3dtke_dyco = .false. !>Convenience pointers integer, pointer :: grid_number @@ -1062,6 +1064,16 @@ module fv_arrays_mod logical :: BCfile_ne_is_open=.false. logical :: BCfile_sw_is_open=.false. end type fv_nest_type + + !3D-SA-TKE (kyf) (modify for data structure) + type sa3dtke_type + real, _ALLOCATABLE :: deform_1(:,:,:) _NULL !< horizontal deformation + real, _ALLOCATABLE :: deform_2(:,:,:) _NULL !< vertical deformation + real, _ALLOCATABLE :: deform_3(:,:,:) _NULL !< 3D TKE transport & pressure correlation + real, _ALLOCATABLE :: dku3d_h(:,:,:) _NULL !< 3D Horizontal Eddy Diffusivity for Momentum + real, _ALLOCATABLE :: dku3d_e(:,:,:) _NULL !< 3D Eddy Diffusivity for TKE + end type sa3dtke_type + !3D-SA-TKE-end type inline_mp_type real, _ALLOCATABLE :: prer(:,:) _NULL @@ -1286,7 +1298,6 @@ module fv_arrays_mod real, _ALLOCATABLE :: va(:,:,:) _NULL real, _ALLOCATABLE :: uc(:,:,:) _NULL !< (uc, vc) are mostly used as the C grid winds real, _ALLOCATABLE :: vc(:,:,:) _NULL - real, _ALLOCATABLE :: ak(:) _NULL real, _ALLOCATABLE :: bk(:) _NULL @@ -1373,6 +1384,7 @@ module fv_arrays_mod integer :: atmos_axes(4) + type(sa3dtke_type) :: sa3dtke_var ! SA-3D TKE (kyf) (modify for data structure) type(inline_mp_type) :: inline_mp type(phys_diag_type) :: phys_diag type(nudge_diag_type) :: nudge_diag @@ -1519,6 +1531,17 @@ subroutine allocate_fv_atmos_type(Atm, isd_in, ied_in, jsd_in, jed_in, is_in, ie allocate ( Atm%va(isd:ied ,jsd:jed ,npz) ) allocate ( Atm%uc(isd:ied+1,jsd:jed ,npz) ) allocate ( Atm%vc(isd:ied ,jsd:jed+1,npz) ) + + !3D-SA-TKE (kyf) (modify for data structure) + if ( Atm%flagstruct%sa3dtke_dyco ) then + allocate ( Atm%sa3dtke_var%deform_1(isd:ied ,jsd:jed ,npz) ) + allocate ( Atm%sa3dtke_var%deform_2(isd:ied ,jsd:jed ,npz) ) + allocate ( Atm%sa3dtke_var%deform_3(isd:ied ,jsd:jed ,npz) ) + allocate ( Atm%sa3dtke_var%dku3d_h(isd:ied ,jsd:jed, npz) ) + allocate ( Atm%sa3dtke_var%dku3d_e(isd:ied ,jsd:jed, npz) ) + endif + !3D-SA-TKE-end + ! For tracer transport: allocate ( Atm%mfx(is:ie+1, js:je, npz) ) allocate ( Atm%mfy(is:ie , js:je+1,npz) ) @@ -1571,6 +1594,16 @@ subroutine allocate_fv_atmos_type(Atm, isd_in, ied_in, jsd_in, jed_in, is_in, ie Atm%va(i,j,k) = real_big Atm%pt(i,j,k) = real_big Atm%delp(i,j,k) = real_big + + !3D-SA-TKE (kyf) (modify for data structure) + if ( Atm%flagstruct%sa3dtke_dyco ) then + Atm%sa3dtke_var%deform_1(i,j,k) = 0. + Atm%sa3dtke_var%deform_2(i,j,k) = 0. + Atm%sa3dtke_var%deform_3(i,j,k) = 0. + Atm%sa3dtke_var%dku3d_h(i,j,k) = 0. + Atm%sa3dtke_var%dku3d_e(i,j,k) = 0. + endif + !3D-SA-TKE-end #ifdef USE_COND Atm%q_con(i,j,k) = 0. #endif diff --git a/model/fv_control.F90 b/model/fv_control.F90 index 921c2ebd0..50d60157f 100644 --- a/model/fv_control.F90 +++ b/model/fv_control.F90 @@ -395,6 +395,9 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split, integer, pointer :: ndims +!The following is used for SA-3D-TKE + logical , pointer :: sa3dtke_dyco + real(kind=R_GRID), pointer :: dx_const real(kind=R_GRID), pointer :: dy_const real(kind=R_GRID), pointer :: deglon_start, deglon_stop, & ! boundaries of latlon patch @@ -954,6 +957,10 @@ subroutine set_namelist_pointers(Atm) external_eta => Atm%flagstruct%external_eta read_increment => Atm%flagstruct%read_increment increment_file_on_native_grid => Atm%flagstruct%increment_file_on_native_grid + +!The following is used for SA-3D-TKE + sa3dtke_dyco => Atm%flagstruct%sa3dtke_dyco + hydrostatic => Atm%flagstruct%hydrostatic phys_hydrostatic => Atm%flagstruct%phys_hydrostatic use_hydro_pressure => Atm%flagstruct%use_hydro_pressure @@ -1083,6 +1090,7 @@ subroutine read_namelist_fv_core_nml(Atm) consv_te, fill, filter_phys, fill_dp, fill_wz, fill_gfs, consv_am, RF_fast, & range_warn, dwind_2d, inline_q, z_tracer, reproduce_sum, adiabatic, do_vort_damp, no_dycore, & tau, tau_w, fast_tau_w_sec, tau_h2o, rf_cutoff, rf_cutoff_w, nf_omega, hydrostatic, fv_sg_adj, sg_cutoff, breed_vortex_inline, & + sa3dtke_dyco, & na_init, nudge_dz, hybrid_z, Make_NH, n_zs_filter, nord_zs_filter, full_zs_filter, reset_eta, & pnats, dnats, dnrts, a2b_ord, remap_t, p_ref, d2_bg_k1, d2_bg_k2, & c2l_ord, dx_const, dy_const, umax, deglat, & diff --git a/model/fv_dynamics.F90 b/model/fv_dynamics.F90 index 36b048136..063589282 100644 --- a/model/fv_dynamics.F90 +++ b/model/fv_dynamics.F90 @@ -156,6 +156,7 @@ module fv_dynamics_mod use fv_regional_mod, only: current_time_in_seconds use boundary_mod, only: nested_grid_BC_apply_intT use fv_arrays_mod, only: fv_grid_type, fv_flags_type, fv_atmos_type, fv_nest_type, fv_diag_type, fv_grid_bounds_type, inline_mp_type + use fv_arrays_mod, only: sa3dtke_type ! for SA-3D-TKE (kyf) (modify for data structure) use fv_nwp_nudge_mod, only: do_adiabatic_init use time_manager_mod, only: get_time @@ -187,8 +188,11 @@ module fv_dynamics_mod subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, & reproduce_sum, kappa, cp_air, zvir, ptop, ks, ncnst, n_split, & - q_split, u, v, w, delz, hydrostatic, pt, delp, q, & + q_split, u, v, w, delz, hydrostatic, & + pt, delp, q, & ps, pe, pk, peln, pkz, phis, q_con, omga, ua, va, uc, vc, & +!The following variable is for SA-3D-TKE (kyf) (modify for data structure) + sa3dtke_var, & ak, bk, mfx, mfy, cx, cy, ze0, hybrid_z, & gridstruct, flagstruct, neststruct, idiag, bd, & parent_grid, domain, diss_est, inline_mp) @@ -255,8 +259,11 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, real, intent(inout) :: vc(bd%isd:bd%ied ,bd%jsd:bd%jed+1,npz) real, intent(inout), dimension(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz):: ua, va - real, intent(in), dimension(npz+1):: ak, bk + !The following variable is for SA-3D-TKE (kyf) (modify for data structure) + real, intent(in), dimension(npz+1):: ak, bk + ! For SA-3D-TKE (kyf) (modify for data structure) + type(sa3dtke_type), intent(inout) :: sa3dtke_var type(inline_mp_type), intent(inout) :: inline_mp ! Accumulated Mass flux arrays: the "Flux Capacitor" @@ -666,7 +673,10 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, #endif grav, hydrostatic, & u, v, w, delz, pt, q, delp, pe, pk, phis, ws, omga, ptop, pfull, ua, va, & - uc, vc, mfx, mfy, cx, cy, pkz, peln, q_con, ak, bk, ks, & + uc, vc, & +!The following variable is for SA-3D-TKE (kyf) (modify for data structure) + sa3dtke_var, & + mfx, mfy, cx, cy, pkz, peln, q_con, ak, bk, ks, & gridstruct, flagstruct, neststruct, idiag, bd, & domain, n_map==1, i_pack, GFDL_interstitial%last_step, diss_est,time_total) call timing_off('DYN_CORE') diff --git a/model/nh_utils.F90 b/model/nh_utils.F90 index 5f3d7c7fd..d2895f168 100644 --- a/model/nh_utils.F90 +++ b/model/nh_utils.F90 @@ -70,6 +70,9 @@ module nh_utils_mod public sim_solver, sim1_solver, sim3_solver public sim3p0_solver, rim_2d public Riem_Solver_c +!This is for sa3dtke + public edge_profile1 + real, parameter:: r3 = 1./3. @@ -1947,6 +1950,73 @@ subroutine edge_profile(q1, q2, q1e, q2e, i1, i2, j1, j2, j, km, dp0, uniform_gr end subroutine edge_profile + subroutine edge_profile1(q1, q1e, i1, i2, km, dp0, limiter) +! Edge profiles for a single scalar quantity + integer, intent(in):: i1, i2 + integer, intent(in):: km + integer, intent(in):: limiter + real, intent(in):: dp0(km) + real, intent(in), dimension(i1:i2,km):: q1 + real, intent(out), dimension(i1:i2,km+1):: q1e +!----------------------------------------------------------------------- + real, dimension(i1:i2,km+1):: qe1, gam ! edge values + real gak(km) + real bet, r2o3, r4o3 + real g0, gk, xt1, xt2, a_bot + integer i, k + +! Assuming grid varying in vertical only + g0 = dp0(2) / dp0(1) + xt1 = 2.*g0*(g0+1. ) + bet = g0*(g0+0.5) + do i=i1,i2 + qe1(i,1) = ( xt1*q1(i,1) + q1(i,2) ) / bet + gam(i,1) = ( 1. + g0*(g0+1.5) ) / bet + enddo + + do k=2,km + gk = dp0(k-1) / dp0(k) + do i=i1,i2 + bet = 2. + 2.*gk - gam(i,k-1) + qe1(i,k) = ( 3.*(q1(i,k-1)+gk*q1(i,k)) - qe1(i,k-1) ) / bet + gam(i,k) = gk / bet + enddo + enddo + + a_bot = 1. + gk*(gk+1.5) + xt1 = 2.*gk*(gk+1.) + do i=i1,i2 + xt2 = gk*(gk+0.5) - a_bot*gam(i,km) + qe1(i,km+1) = ( xt1*q1(i,km) + q1(i,km-1) - a_bot*qe1(i,km) ) / xt2 + enddo + + do k=km,1,-1 + do i=i1,i2 + qe1(i,k) = qe1(i,k) - gam(i,k)*qe1(i,k+1) + enddo + enddo + +!------------------ +! Apply constraints +!------------------ + if ( limiter/=0 ) then ! limit the top & bottom winds + do i=i1,i2 +! Top + if ( q1(i,1)*qe1(i,1) < 0. ) qe1(i,1) = 0. +! Surface: + if ( q1(i,km)*qe1(i,km+1) < 0. ) qe1(i,km+1) = 0. + enddo + endif + + do k=1,km+1 + do i=i1,i2 + q1e(i,k) = qe1(i,k) + enddo + enddo + + end subroutine edge_profile1 + + subroutine edge_profile_0grad(q1, q2, q1e, q2e, i1, i2, j1, j2, j, km, dp0, uniform_grid, limiter) ! Optimized for wind profile reconstruction: ! Added this option by Henry Juang and Xiaqiong Zhou 1/21/2021 diff --git a/model/sw_core.F90 b/model/sw_core.F90 index 6ed4a2e50..54d7e7bd7 100644 --- a/model/sw_core.F90 +++ b/model/sw_core.F90 @@ -525,7 +525,9 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, & zvir, sphum, nq, q, k, km, inline_q, & dt, hord_tr, hord_mt, hord_vt, hord_tm, hord_dp, nord, & nord_v, nord_w, nord_t, dddmp, d2_bg, d4_bg, damp_v, damp_w, & - damp_t, d_con, hydrostatic, gridstruct, flagstruct, bd) + damp_t, d_con, hydrostatic, gridstruct, flagstruct, bd, & +!The following optional variable is for SA-3D-TKE (kyf) + dku3d_h) integer, intent(IN):: hord_tr, hord_mt, hord_vt, hord_tm, hord_dp integer, intent(IN):: nord !< nord=1 divergence damping; (del-4) or 3 (del-8) @@ -537,6 +539,7 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, & real , intent(IN):: zvir real, intent(in):: damp_v, damp_w, damp_t, kgb type(fv_grid_bounds_type), intent(IN) :: bd + real, intent(inout):: divg_d(bd%isd:bd%ied+1,bd%jsd:bd%jed+1) !< divergence real, intent(IN), dimension(bd%isd:bd%ied, bd%jsd:bd%jed):: z_rat real, intent(INOUT), dimension(bd%isd:bd%ied, bd%jsd:bd%jed):: delp, pt, ua, va @@ -559,6 +562,8 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, & real, intent(OUT), dimension(bd%isd:bd%ied,bd%js:bd%je+1):: cry_adv, yfx_adv type(fv_grid_type), intent(IN), target :: gridstruct type(fv_flags_type), intent(IN), target :: flagstruct +!The following optional variable is for SA-3D-TKE (kyf) + real, intent(IN), optional, dimension(bd%isd:bd%ied,bd%jsd:bd%jed):: dku3d_h ! Local: logical:: sw_corner, se_corner, ne_corner, nw_corner real :: ut(bd%isd:bd%ied+1,bd%jsd:bd%jed) @@ -1453,7 +1458,13 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, & do j=js,je+1 do i=is,ie+1 +!The following is for SA-3D-TKE + if(flagstruct%sa3dtke_dyco) then + damp2 = abs(dt)*dku3d_h(i,j) + else damp2 = gridstruct%da_min_c*max(d2_bg, min(0.20, dddmp*vort(i,j))) ! del-2 + endif + vort(i,j) = damp2*delpc(i,j) + dd8*divg_d(i,j) ke(i,j) = ke(i,j) + vort(i,j) enddo