Skip to content

Commit

Permalink
more on diagn
Browse files Browse the repository at this point in the history
  • Loading branch information
@oksanaguba
oksanaguba committed Dec 17, 2024
1 parent 1e20422 commit bcb8ce4
Show file tree
Hide file tree
Showing 3 changed files with 113 additions and 25 deletions.
5 changes: 5 additions & 0 deletions components/homme/src/theta-l/element_state.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -89,6 +89,11 @@ module element_state

#ifdef ENERGY_DIAGNOSTICS
! Energy equation:
real (kind=real_kind) :: PE(np,np)
real (kind=real_kind) :: PEexpected(np,np)
real (kind=real_kind) :: ieterm1(np,np)
real (kind=real_kind) :: keterm1(np,np)

real (kind=real_kind) :: KEu_horiz1(np,np)
real (kind=real_kind) :: KEu_horiz2(np,np)
real (kind=real_kind) :: KEu_vert1(np,np)
Expand Down
102 changes: 77 additions & 25 deletions components/homme/src/theta-l/share/prim_advance_mod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1149,7 +1149,9 @@ subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,&
real (kind=real_kind), dimension(np,np,nlev) :: vgrad_p

real (kind=real_kind) :: temp(np,np,nlev)
real (kind=real_kind) :: tempp(np,np,nlevp)
real (kind=real_kind) :: vtemp(np,np,2,nlev) ! generic gradient storage
real (kind=real_kind) :: w_m(np,np,nlev)
real (kind=real_kind), dimension(np,np) :: sdot_sum ! temporary field
real (kind=real_kind) :: v1,v2,w,d_eta_dot_dpdn_dn, T0
integer :: i,j,k,kptr,ie, nlyr_tot
Expand Down Expand Up @@ -1181,6 +1183,8 @@ subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,&
vtheta(:,:,:) = vtheta_dp(:,:,:)/dp3d(:,:,:)
phi_i => elem(ie)%state%phinh_i(:,:,:,n0)

!make w_i alias too

#ifdef HOMMEDA
!repeated code
rheighti = phi_i/g + r0
Expand Down Expand Up @@ -1485,6 +1489,10 @@ subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,&
phi_tens(:,:,k) = (-phi_vadv_i(:,:,k) - v_gradphinh_i(:,:,k))*scale1 &
+ scale1*g*elem(ie)%state%w_i(:,:,k,n0)

#if defined HOMMEDA && defined ENERGY_DIAGNOSTICS
phi_tens_notopo(:,:,k) = phi_tens(:,:,k)
#endif


! ================================================
! v1,v2 tendencies:
Expand All @@ -1496,7 +1504,21 @@ subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,&
call i2m(elem(ie)%state%w_i(:,:,:,n0)*gradw_i(:,:,2,:),temp(:,:,:))
vtemp(:,:,2,:) = temp(:,:,:)

call i2m(dpnh_dp_i*gradphinh_i(:,:,1,:),temp(:,:,:))
#ifdef HOMMEDA
mgrad(:,:,1,:) = temp(:,:,:)*invrhatm(:,:,:)
#else
mgrad(:,:,1,:) = temp(:,:,:)
#endif
call i2m(dpnh_dp_i*gradphinh_i(:,:,2,:),temp(:,:,:))
#ifdef HOMMEDA
mgrad(:,:,2,:) = temp(:,:,:)*invrhatm(:,:,:)
#else
mgrad(:,:,2,:) = temp(:,:,:)
#endif

call i2m(elem(ie)%state%w_i(:,:,:,n0)*elem(ie)%state%w_i(:,:,:,n0),temp)
call i2m(elem(ie)%state%w_i(:,:,:,n0),w_m)

do k=1,nlev
! theta - tendency on levels
Expand Down Expand Up @@ -1534,18 +1556,12 @@ subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,&
#ifdef HOMMEDA
wvor(:,:,1,k) = wvor(:,:,1,k) * invrhatm(:,:,k)
wvor(:,:,2,k) = wvor(:,:,2,k) * invrhatm(:,:,k)
#endif

wvor(:,:,1,k) = wvor(:,:,1,k) - vtemp(:,:,1,k) * invrhatm(:,:,k)
wvor(:,:,2,k) = wvor(:,:,2,k) - vtemp(:,:,2,k) * invrhatm(:,:,k)
#else
wvor(:,:,1,k) = wvor(:,:,1,k) - vtemp(:,:,1,k)
wvor(:,:,2,k) = wvor(:,:,2,k) - vtemp(:,:,2,k)

!there is already a DA correction in gradw_i
! wvor(:,:,1,k) = wvor(:,:,1,k) - (elem(ie)%state%w_i(:,:,k,n0)*gradw_i(:,:,1,k) +&
! elem(ie)%state%w_i(:,:,k+1,n0)*gradw_i(:,:,1,k+1))/2
! wvor(:,:,2,k) = wvor(:,:,2,k) - (elem(ie)%state%w_i(:,:,k,n0)*gradw_i(:,:,2,k) +&
! elem(ie)%state%w_i(:,:,k+1,n0)*gradw_i(:,:,2,k+1))/2


#endif

KE(:,:,k) = ( elem(ie)%state%v(:,:,1,k,n0)**2 + elem(ie)%state%v(:,:,2,k,n0)**2)/2
gradKE(:,:,:,k) = gradient_sphere(KE(:,:,k),deriv,elem(ie)%Dinv)
Expand Down Expand Up @@ -1588,14 +1604,6 @@ subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,&
gradexner(:,:,2,k) = gradexner(:,:,2,k)*(Rgas/Cp)*exner(:,:,k)/pnh(:,:,k)
#endif

!gradphinh_i already has DA correction
! special averaging of dpnh/dpi grad(phi) for E conservation
mgrad(:,:,1,k) = (dpnh_dp_i(:,:,k)*gradphinh_i(:,:,1,k)+ &
dpnh_dp_i(:,:,k+1)*gradphinh_i(:,:,1,k+1))/2
mgrad(:,:,2,k) = (dpnh_dp_i(:,:,k)*gradphinh_i(:,:,2,k)+ &
dpnh_dp_i(:,:,k+1)*gradphinh_i(:,:,2,k+1))/2


!OG do pgrad DA later !
if (pgrad_correction==1) then
T0 = TREF-tref_lapse_rate*TREF*Cp/g ! = 97
Expand Down Expand Up @@ -1647,9 +1655,8 @@ subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,&
-wvor(i,j,2,k) )*scale1

#ifdef HOMMEDA
vtens1(i,j,k) = vtens1(i,j,k) - scale1*(elem(ie)%state%w_i(i,j,k,n0)+elem(ie)%state%w_i(i,j,k+1,n0))/2*( v1/rheightm(i,j,k) &
+ elem(ie)%fcorcosine(i,j) )
vtens2(i,j,k) = vtens2(i,j,k) - scale1*(elem(ie)%state%w_i(i,j,k,n0)+elem(ie)%state%w_i(i,j,k+1,n0))/2*v2/rheightm(i,j,k)
vtens1(i,j,k) = vtens1(i,j,k) - scale1*w_m(i,j,k)*( v1/rheightm(i,j,k) + elem(ie)%fcorcosine(i,j) )
vtens2(i,j,k) = vtens2(i,j,k) - scale1*w_m(i,j,k)*v2/rheightm(i,j,k)
#endif

#endif
Expand All @@ -1665,7 +1672,12 @@ subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,&
! diagnostics. not performance critical, dont thread
! =========================================================
if (compute_diagnostics) then

elem(ie)%accum%PE=0
elem(ie)%accum%PEexpected=0

elem(ie)%accum%ieterm1=0
elem(ie)%accum%keterm1=0

elem(ie)%accum%KEu_horiz1=0
elem(ie)%accum%KEu_horiz2=0
Expand All @@ -1690,6 +1702,39 @@ subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,&
elem(ie)%accum%P1=0
elem(ie)%accum%P2=0

!use temp and vtemp
call m2i(divdp,tempp)

!divdp = -dp3d_tens
do k=2,nlev
elem(ie)%accum%PE(:,:) = elem(ie)%accum%PE(:,:) + phi_i(:,:,k)*tempp(:,:,k) &
- dp3d_i(:,:,k)*phi_tens_notopo(:,:,k)
elem(ie)%accum%PEexpected(:,:) = elem(ie)%accum%PEexpected(:,:) + g*dp3d_i(:,:,k)*elem(ie)%state%w_i(:,:,k,n0)
enddo
do k=1,nlevp,nlev
elem(ie)%accum%PE(:,:) = elem(ie)%accum%PE(:,:) + phi_i(:,:,k)*tempp(:,:,k)/2 &
- dp3d_i(:,:,k)*phi_tens_notopo(:,:,k)/2
elem(ie)%accum%PEexpected(:,:) = elem(ie)%accum%PEexpected(:,:) + g*dp3d_i(:,:,k)*elem(ie)%state%w_i(:,:,k,n0)/2
enddo

!note that mgrad is overwritten above, so do not re-use, recompute
call i2m(dpnh_dp_i*gradphinh_i(:,:,1,:),temp(:,:,:))
mgrad(:,:,1,:) = temp(:,:,:)
call i2m(dpnh_dp_i*gradphinh_i(:,:,2,:),temp(:,:,:))
mgrad(:,:,2,:) = temp(:,:,:)
do k=2,nlev
elem(ie)%accum%ieterm1(:,:) = elem(ie)%accum%ieterm1(:,:) + dpnh_dp_i(:,:,k)*dp3d_i(:,:,k)*v_gradphinh_i(:,:,k)
enddo
do k=1,nlevp,nlev
elem(ie)%accum%ieterm1(:,:) = elem(ie)%accum%ieterm1(:,:) + dpnh_dp_i(:,:,k)*dp3d_i(:,:,k)*v_gradphinh_i(:,:,k)/2
enddo
do k=1,nlev
elem(ie)%accum%keterm1(:,:) = elem(ie)%accum%keterm1(:,:) - dp3d(:,:,k)*invrhatm(:,:,k)* &
( elem(ie)%state%v(:,:,1,k,n0)*mgrad(:,:,1,k) + elem(ie)%state%v(:,:,2,k,n0)*mgrad(:,:,2,k) )
enddo



do k =1,nlev
do j=1,np
do i=1,np
Expand Down Expand Up @@ -1762,12 +1807,19 @@ subroutine compute_andor_apply_rhs(np1,nm1,n0,dt2,elem,hvcoord,hybrid,&

!done
! Form PEhoriz2
! elem(ie)%accum%PEhoriz2(i,j)=elem(ie)%accum%PEhoriz2(i,j) &
! -dp3d(i,j,k)* &
! (elem(ie)%state%v(i,j,1,k,n0)* &
! (gradphinh_i(i,j,1,k)+gradphinh_i(i,j,1,k+1))/2 + &
! elem(ie)%state%v(i,j,2,k,n0)* &
! (gradphinh_i(i,j,2,k)+gradphinh_i(i,j,2,k+1))/2 )

elem(ie)%accum%PEhoriz2(i,j)=elem(ie)%accum%PEhoriz2(i,j) &
-dp3d(i,j,k)* &
(elem(ie)%state%v(i,j,1,k,n0)* &
(gradphinh_i(i,j,1,k)+gradphinh_i(i,j,1,k+1))/2 + &
elem(ie)%state%v(i,j,2,k,n0)* &
(gradphinh_i(i,j,2,k)+gradphinh_i(i,j,2,k+1))/2 )
(v_gradphinh_i(i,j,k) + v_gradphinh_i(i,j,k+1))/2




!done as these are 0
! Form PEvert1
Expand Down
31 changes: 31 additions & 0 deletions components/homme/src/theta-l/share/prim_state_mod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -148,6 +148,9 @@ subroutine prim_printstate(elem, tl,hybrid,hvcoord,nets,nete)
real (kind=real_kind) :: PEhorz1,PEhorz2
real (kind=real_kind) :: KEH1,KEH2,KEV1,KEV2
real (kind=real_kind) :: KEwH1,KEwH2,KEwH3,KEwV1,KEwV2

real (kind=real_kind) :: PEscalar, PEexpected_scalar, iet1s, ket1s

real (kind=real_kind) :: ddt_tot,ddt_diss, ddt_diss_adj
integer :: n0, n0q
integer :: npts,n,q
Expand All @@ -169,6 +172,9 @@ subroutine prim_printstate(elem, tl,hybrid,hvcoord,nets,nete)
PEner = 0
IEner = 0
muvalue = 0

PEscalar = 0; PEexpected_scalar = 0; iet1s = 0; ket1s = 0;

! dynamics timelevels
n0=tl%n0
call TimeLevel_Qdp(tl, qsplit, n0q) ! get n0 level into n0q
Expand Down Expand Up @@ -638,6 +644,24 @@ subroutine prim_printstate(elem, tl,hybrid,hvcoord,nets,nete)
PEvert2 = global_integral(elem, tmp(:,:,nets:nete),hybrid,npts,nets,nete)
PEvert2 = PEvert2*scale

do ie=nets,nete
tmp(:,:,ie) = elem(ie)%accum%PE
enddo
PEscalar = global_integral(elem, tmp(:,:,nets:nete),hybrid,npts,nets,nete)*scale
do ie=nets,nete
tmp(:,:,ie) = elem(ie)%accum%PEexpected
enddo
PEexpected_scalar = global_integral(elem, tmp(:,:,nets:nete),hybrid,npts,nets,nete)*scale

do ie=nets,nete
tmp(:,:,ie) = elem(ie)%accum%ieterm1
enddo
iet1s = global_integral(elem, tmp(:,:,nets:nete),hybrid,npts,nets,nete)*scale
do ie=nets,nete
tmp(:,:,ie) = elem(ie)%accum%keterm1
enddo
ket1s = global_integral(elem, tmp(:,:,nets:nete),hybrid,npts,nets,nete)*scale

! KE->IE
do ie=nets,nete
tmp(:,:,ie) = elem(ie)%accum%T01
Expand Down Expand Up @@ -711,6 +735,7 @@ subroutine prim_printstate(elem, tl,hybrid,hvcoord,nets,nete)
write(iulog,'(a,2e22.14)')'KEu v-adv,sum=0:',KEV1,KEV2
write(iulog,'(a,3e22.14)')'IE v-adv,sum=0:',IEvert1,PEvert1
write(iulog,'(a,2e22.14)')'KE->I+P,I+P->KE:',(T1+PEhorz2),(S1+PEhorz1)
write(iulog,'(a,2e22.14)')'PEhorz1,PEhorz2:',PEhorz1,PEhorz2

ddt_tot = (KEner(2)-KEner(1))/dt
ddt_diss = ddt_tot -(T1+PEhorz2)
Expand All @@ -731,6 +756,12 @@ subroutine prim_printstate(elem, tl,hybrid,hvcoord,nets,nete)
write(iulog,'(a,3e22.14)')'KEw h-adv,sum=0:',KEwH1+KEwH3,KEwH2
write(iulog,'(a,2e22.14)')'KEu v-adv,sum=0:',KEV1,KEV2
write(iulog,'(a,2e22.14)')'KEw v-adv,sum=0:',KEwV1,KEwV2

write(iulog,'(a,2e22.14)')'PEhorz1,PEhorz2:',PEhorz1,PEhorz2

write(iulog,'(a,2e22.14)')'PE_t,PE_t exp:',PEscalar, PEexpected_scalar
write(iulog,'(a,2e22.14)')'ieterm1,keterm1:', iet1s, ket1s

write(iulog,'(a,2e22.14)')'PE h-adv, sum=0:',PEhorz1,PEhorz2
write(iulog,'(a,2e22.14)')'PE v-adv, sum=0:',PEvert1,PEvert2
write(iulog,'(a,3e22.14)')'IE v-adv, sum=0:',IEvert1,IEvert2
Expand Down

0 comments on commit bcb8ce4

Please sign in to comment.