Skip to content
Merged
Show file tree
Hide file tree
Changes from 2 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 @@ -2195,9 +2195,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 @@ -1915,6 +1915,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