Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 3 additions & 3 deletions driver/UFS/atmosphere.F90
Original file line number Diff line number Diff line change
Expand Up @@ -2197,9 +2197,9 @@ subroutine atmos_phys_driver_statein (IPD_Control, IPD_Statein, IPD_Tbd, Atm_blo
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)))
IPD_Statein%def_1(im,k) = _DBL_(_RL_(Atm(mygrid)%sa3dtke_var%deform_1(i,j,k1)))
IPD_Statein%def_2(im,k) = _DBL_(_RL_(Atm(mygrid)%sa3dtke_var%deform_2(i,j,k1)))
IPD_Statein%def_3(im,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)))
Expand Down
158 changes: 75 additions & 83 deletions model/dyn_core.F90
Original file line number Diff line number Diff line change
Expand Up @@ -897,7 +897,6 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp,
enddo
enddo
endif

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), &
Expand Down Expand Up @@ -2206,25 +2205,16 @@ subroutine diff3d(npx, npy, npz, nq, ua, va, w, q, &
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)
real:: ut(bd%isd:bd%ied,bd%jsd:bd%jed)
real:: vt(bd%isd:bd%ied,bd%jsd:bd%jed)

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 :: tke_1(bd%isd:bd%ied,bd%jsd:bd%jed)
real :: kqdedy(bd%isd:bd%ied,bd%jsd:bd%jed,npz)
real :: dkqdedydy(bd%isd:bd%ied,bd%jsd:bd%jed,npz)
real :: kqdedx(bd%isd:bd%ied,bd%jsd:bd%jed,npz)
real :: dkqdedxdx(bd%isd:bd%ied,bd%jsd:bd%jed,npz)
real :: tmp, tmp1, tmp2

real :: dz

Expand Down Expand Up @@ -2252,61 +2242,61 @@ subroutine diff3d(npx, npy, npz, nq, ua, va, w, q, &
! 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 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 j=js-1,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
do i=is-1,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))
dudy(i,j,k) = 0.5*rarea(i,j)*(vt(i,j+1)-vt(i,j-1))
dvdx(i,j,k) = 0.5*rarea(i,j)*(ut(i+1,j)-ut(i-1,j))
enddo
enddo
!-------------------------------------
! get du/dx and dv/dy
do j=js,je
do i=is,ie+1
do i=is-1,ie+1
ut(i,j) = ua(i,j,k)*dy(i,j)
enddo
enddo
do j=js,je+1
do j=js-1,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))
dudx(i,j,k) = 0.5*rarea(i,j)*(ut(i+1,j)-ut(i-1,j))
dvdy(i,j,k) = 0.5*rarea(i,j)*(vt(i,j+1)-vt(i,j-1))
enddo
enddo
!-------------------------------------
! get dw/dx and dw/dy
do j=js,je
do i=is,ie+1
do i=is-1,ie+1
ut(i,j) = w(i,j,k)*dy(i,j)
enddo
enddo
do j=js,je+1
do j=js-1,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))
dwdx(i,j,k) = 0.5*rarea(i,j)*(ut(i+1,j)-ut(i-1,j))
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please check closely where dwdx, etc. are computed and whether a full-grid-length centered difference is needed here. This will compute dwdx at cell centers instead of cell interfaces, and necessitates the additional interpolation of dwdx below on lines 2338--2357.

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please check closely where dwdx, etc. are computed and whether a full-grid-length centered difference is needed here. This will compute dwdx at cell centers instead of cell interfaces, and necessitates the additional interpolation of dwdx below on lines 2338--2357.

@lharris4 yes, dwdx is computed at cell centers (vertically at layers) and then 2 layer vertical average of dwdx is used for deform_1 computation which is defined at the interfaces

dwdy(i,j,k) = 0.5*rarea(i,j)*(vt(i,j+1)-vt(i,j-1))
enddo
enddo
enddo !z loop
Expand All @@ -2316,34 +2306,55 @@ subroutine diff3d(npx, npy, npz, nq, ua, va, w, q, &

! 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
!$OMP zh,dudz,dvdz,dwdz) &
!$OMP private(dz)
do k=2,npz
do j=js,je
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
dz = 0.5*(zh(i,j,k-1) - zh(i,j,k+1))
dudz(i,j,k) = (ua(i,j,k-1)-ua(i,j,k))/dz
dvdz(i,j,k) = (va(i,j,k-1)-va(i,j,k))/dz
dwdz(i,j,k) = (w(i,j,k-1)-w(i,j,k))/dz
enddo
enddo
enddo
do j=js,je
do i=is,ie
dudz(i,j,1) = 0.
dvdz(i,j,1) = 0.
dwdz(i,j,1) = 0.
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)
!$OMP private(tmp,tmp1,tmp2)
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)
if(k > 1) then
tmp1 = 0.5*(dwdx(i,j,k-1)+dwdx(i,j,k))
tmp2 = 0.5*(dwdy(i,j,k-1)+dwdy(i,j,k))
tmp = tmp1*dudz(i,j,k)+tmp2*dvdz(i,j,k)
else
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why are the branches for the top and bottom layers needed?

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why are the branches for the top and bottom layers needed?

@lharris4 Actually, the deform_1 value at k=1 (npz in the CCPP 3D TKE scheme) is not used in the 3D TKE scheme, but we still assign a reasonable deform_1 value at k=1 for a safety.

tmp1 = dwdx(i,j,k)
tmp2 = dwdy(i,j,k)
tmp = tmp1*dudz(i,j,k)+tmp2*dvdz(i,j,k)
endif
deform_1(i,j,k)= 2*dwdz(i,j,k)**2+ &
dudz(i,j,k)**2+dvdz(i,j,k)**2+tmp
if(k < npz) then
tmp1 = 0.5*(dudz(i,j,k)+dudz(i,j,k+1))
tmp2 = 0.5*(dvdz(i,j,k)+dvdz(i,j,k+1))
tmp = tmp1*dwdx(i,j,k)+tmp2*dwdy(i,j,k)
else
tmp1 = dudz(i,j,k)
tmp2 = dvdz(i,j,k)
tmp = tmp1*dwdx(i,j,k)+tmp2*dwdy(i,j,k)
endif
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+ &
Expand All @@ -2363,87 +2374,68 @@ subroutine diff3d(npx, npy, npz, nq, ua, va, w, q, &
! 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)
!$OMP dku3d_e,kqdedy,dkqdedydy,kqdedx,dkqdedxdx,sgs_tke) &
!$OMP private(ut,vt,tke_1,tmp)
do k=1,npz
!-------------------------------------
do j=js,je+2
do i=is,ie+2
do j=js-1,je+1
do i=is-1,ie+1
tke_1(i,j)=q(i,j,k,sgs_tke)
enddo
enddo
do j=js,je+2
do j=js-1,je+1
do i=is,ie
vt(i,j)=tke_1(i,j)*dx(i,j)
enddo
enddo
do j=js,je+1
do j=js-1,je
do i=is,ie
dedy_1(i,j,k)=rarea(i,j)*(vt(i,j+1)-vt(i,j))
tmp = 0.5*(dku3d_e(i,j,k)+dku3d_e(i,j+1,k))
kqdedy(i,j,k)=rarea(i,j)*tmp*(vt(i,j+1)-vt(i,j))
enddo
enddo
do j=js,je+1
do j=js-1,je
do i=is,ie
vt(i,j)=dedy_1(i,j,k)*dx(i,j)
vt(i,j)=kqdedy(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))
dkqdedydy(i,j,k)=rarea(i,j)*(vt(i,j)-vt(i,j-1))
enddo
enddo
!-------------------------------------
do j=js,je
do i=is,ie+2
do i=is-1,ie+1
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))
do i=is-1,ie
tmp = 0.5*(dku3d_e(i,j,k)+dku3d_e(i+1,j,k))
kqdedx(i,j,k)=rarea(i,j)*tmp*(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)
do i=is-1,ie
ut(i,j)=kqdedx(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))
dkqdedxdx(i,j,k)=rarea(i,j)*(ut(i,j)-ut(i-1,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)
!$OMP dkqdedxdx,dkqdedydy)
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)
deform_3(i,j,k)=dkqdedxdx(i,j,k)+dkqdedydy(i,j,k)
enddo
enddo
enddo
Expand Down
9 changes: 9 additions & 0 deletions model/fv_arrays.F90
Original file line number Diff line number Diff line change
Expand Up @@ -1921,6 +1921,15 @@ subroutine deallocate_fv_atmos_type(Atm)
deallocate ( Atm%va )
deallocate ( Atm%uc )
deallocate ( Atm%vc )
!3D-SA-TKE (kyf) (modify for data structure)
if ( Atm%flagstruct%sa3dtke_dyco ) then
deallocate ( Atm%sa3dtke_var%deform_1 )
deallocate ( Atm%sa3dtke_var%deform_2 )
deallocate ( Atm%sa3dtke_var%deform_3 )
deallocate ( Atm%sa3dtke_var%dku3d_h )
deallocate ( Atm%sa3dtke_var%dku3d_e )
endif
!3D-SA-TKE-end
deallocate ( Atm%mfx )
deallocate ( Atm%mfy )
deallocate ( Atm%cx )
Expand Down
5 changes: 2 additions & 3 deletions model/sw_core.F90
Original file line number Diff line number Diff line change
Expand Up @@ -1458,11 +1458,10 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, &

do j=js,je+1
do i=is,ie+1
damp2 = gridstruct%da_min_c*max(d2_bg, min(0.20, dddmp*vort(i,j))) ! del-2
!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
damp2 = max(damp2, abs(dt)*dku3d_h(i,j))
endif

vort(i,j) = damp2*delpc(i,j) + dd8*divg_d(i,j)
Expand Down