Skip to content

changes to owen's branch #13

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 8 commits into
base: og/copy-owen/deep_atm_10_15_23
Choose a base branch
from
Open

changes to owen's branch #13

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
2ee8cca
dry test changes
oksanaguba May 27, 2024
943be52
add back g for compatibility
oksanaguba May 29, 2024
63271a5
revert IC to master fork
oksanaguba May 29, 2024
bcfc2cb
revert is_deep type to match baro mod
oksanaguba May 29, 2024
f104032
reorg eos calls to have both versions
oksanaguba May 29, 2024
79eb06d
two options for imex
oksanaguba May 29, 2024
dedc4e7
comments
oksanaguba May 29, 2024
547e4ce
typo
oksanaguba May 29, 2024
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
1 change: 1 addition & 0 deletions components/homme/src/share/physical_constants.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,7 @@ module physical_constants
real (kind=real_kind), public, parameter :: rearth0 = 6.376D6 ! m
real (kind=real_kind), public :: rearth = rearth0 ! m
real (kind=real_kind), public, parameter :: gravit = 9.80616D0 ! m s^-2
real (kind=real_kind), public, parameter :: g = 9.80616D0 ! m s^-2
real (kind=real_kind), public :: gravitinv = 1.0_real_kind/gravit
real (kind=real_kind), public, parameter :: omega0 = 7.292D-5 ! s^-1
real (kind=real_kind), public :: omega = omega0
Expand Down
22 changes: 11 additions & 11 deletions components/homme/src/test_src/dcmip12_wrapper.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,7 +91,7 @@ subroutine dcmip2012_test1_1(elem,hybrid,hvcoord,nets,nete,time,n0,n1)

dp = pressure_thickness(ps,k,hvcoord)
call set_state(u,v,w,T,ps,phis,p,dp,zm(k),g, i,j,k,elem(ie),n0,n1)
if(time==0) call set_tracers(q,qsize,dp,zm(k),i,j,k,lat,lon,elem(ie))
if(time==0) call set_tracers(q,qsize,dp,i,j,k,lat,lon,elem(ie))

enddo; enddo; enddo; enddo

Expand Down Expand Up @@ -167,7 +167,7 @@ subroutine dcmip2012_test1_1_conv(elem,hybrid,hvcoord,nets,nete,time,n0,n1)

dp = pressure_thickness(ps,k,hvcoord)
call set_state(u,v,w,T,ps,phis,p,dp,zm(k),g, i,j,k,elem(ie),n0,n1)
if(time==0) call set_tracers(q,qsize,dp,zm(k),i,j,k,lat,lon,elem(ie))
if(time==0) call set_tracers(q,qsize,dp,i,j,k,lat,lon,elem(ie))

enddo; enddo; enddo; enddo

Expand Down Expand Up @@ -235,7 +235,7 @@ subroutine dcmip2012_test1_2(elem,hybrid,hvcoord,nets,nete,time,n0,n1)
call test1_advection_hadley(time,lon,lat,p,z,zcoords,u,v,w,t,phis,ps,rho,q(1),q(2))
dp = pressure_thickness(ps,k,hvcoord)
call set_state(u,v,w,T,ps,phis,p,dp,zm(k),g, i,j,k,elem(ie),n0,n1)
if(time==0) call set_tracers(q,qsize,dp,zm(k),i,j,k,lat,lon,elem(ie))
if(time==0) call set_tracers(q,qsize,dp,i,j,k,lat,lon,elem(ie))

enddo; enddo; enddo; enddo

Expand Down Expand Up @@ -306,7 +306,7 @@ subroutine dcmip2012_test1_3(elem,hybrid,hvcoord,nets,nete,time,n0,n1,deriv)
call test1_advection_orography(lon,lat,p,z,zcoords,cfv,use_eta,hyam,hybm,gc,u,v,w,t,phis,ps,rho,q(1),q(2),q(3),q(4))
dp = pressure_thickness(ps,k,hvcoord)
call set_state(u,v,w,T,ps,phis,p,dp,zm(k),g, i,j,k,elem(ie),n0,n1)
if(time==0) call set_tracers(q,qsize,dp,zm(k),i,j,k,lat,lon,elem(ie))
if(time==0) call set_tracers(q,qsize,dp,i,j,k,lat,lon,elem(ie))

enddo; enddo; enddo; enddo

Expand Down Expand Up @@ -373,7 +373,7 @@ subroutine dcmip2012_test2_0(elem,hybrid,hvcoord,nets,nete)
!let's get an analytical \phi
he = (T0 - T)/gamma
call set_state(u,v,w,T,ps,phis,p,dp,he,g, i,j,k,elem(ie),1,nt)
call set_tracers(q,qsize,dp,he,i,j,k,lat,lon,elem(ie))
call set_tracers(q,qsize,dp,i,j,k,lat,lon,elem(ie))
enddo; enddo; enddo;
do k=1,nlevp; do j=1,np; do i=1,np
call get_coordinates(lat,lon,hyai,hybi, i,j,k,elem(ie),hvcoord)
Expand Down Expand Up @@ -432,7 +432,7 @@ subroutine dcmip2012_test2_x(elem,hybrid,hvcoord,nets,nete,shear)
! call set_state(u,v,w,T,ps,phis,p,dp,zm(k),g, i,j,k,elem(ie),1,nt)
! This test obtains analytical height and returns it, so, we use it for \phi ...
call set_state(u,v,w,T,ps,phis,p,dp,z,g, i,j,k,elem(ie),1,nt)
call set_tracers(q,qsize,dp,z,i,j,k,lat,lon,elem(ie))
call set_tracers(q,qsize,dp,i,j,k,lat,lon,elem(ie))
! ... or we can use discrete hydro state to init \phi.

enddo; enddo; enddo;
Expand Down Expand Up @@ -602,7 +602,7 @@ subroutine dcmip2012_test3(elem,hybrid,hvcoord,nets,nete)
call test3_gravity_wave(lon,lat,p,z,zcoords,use_eta,hyam,hybm,u,v,w,T,T_mean,phis,ps,rho,rho_mean,q(1))
dp = pressure_thickness(ps,k,hvcoord)
call set_state(u,v,w,T,ps,phis,p,dp,zm(k),g, i,j,k,elem(ie),1,nt)
call set_tracers(q,qsize, dp,z,i,j,k,lat,lon,elem(ie))
call set_tracers(q,qsize, dp,i,j,k,lat,lon,elem(ie))
enddo; enddo; enddo;
do k=1,nlevp; do j=1,np; do i=1,np
call get_coordinates(lat,lon,hyai,hybi, i,j,k,elem(ie),hvcoord)
Expand Down Expand Up @@ -662,7 +662,7 @@ subroutine dcmip2012_test4_init(elem,hybrid,hvcoord,nets,nete)

!init only <=qsize tracers
qs = min(qsize,3)
call set_tracers(qarray(1:qs),qs, dp,z,i,j,k,lat,lon,elem(ie))
call set_tracers(qarray(1:qs),qs, dp,i,j,k,lat,lon,elem(ie))
enddo; enddo; enddo;

do k=1,nlevp
Expand Down Expand Up @@ -785,19 +785,19 @@ real(rl) function pressure_thickness(ps,k,hv)


!_____________________________________________________________________
subroutine set_tracers(q,nq, dp,zm,i,j,k,lat,lon,elem)
subroutine set_tracers(q,nq, dp,i,j,k,lat,lon,elem)
use physical_constants, only: rearth
use deep_atm_mod, only: r_hat_from_phi, phi_from_z
! set tracer values at node(i,j,k)

real(rl), intent(in) :: q(nq), dp,zm, lat, lon
real(rl), intent(in) :: q(nq), dp, lat, lon
integer, intent(in) :: i,j,k,nq
type(element_t), intent(inout) :: elem
real(rl) :: r_hat, pdensity
real(rl), parameter :: wl = 1.0 ! checkerboard wavelength in dg
integer :: qi

r_hat = r_hat_from_phi(phi_from_z(zm))
!r_hat = r_hat_from_phi(phi_from_z(zm))
pdensity = dp

if (nq>qsize) call abortmp('qsize set too small for dcmip test case')
Expand Down
181 changes: 29 additions & 152 deletions components/homme/src/test_src/dcmip16_wrapper.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,8 @@ subroutine dcmip2016_test1(elem,hybrid,hvcoord,nets,nete)
integer, intent(in) :: nets,nete ! start, end element index

integer, parameter :: use_zcoords = 0 ! use vertical pressure coordinates
logical, parameter :: is_deep = atm_is_deep ! use shallow atmosphere approximation
!logical, parameter :: is_deep = atm_is_deep ! use shallow atmosphere approximation
integer, parameter :: is_deep = 1 ! use shallow atmosphere approximation
integer, parameter :: pertt = 0 ! use exponential perturbation type
real(rl), parameter :: dcmip_X = 1.0_rl ! full scale planet
integer :: moist ! use moist version
Expand Down Expand Up @@ -115,172 +116,48 @@ subroutine dcmip2016_test1(elem,hybrid,hvcoord,nets,nete)

if (qsize<6) call abortmp('ERROR: test requires qsize>=6')

! set initial conditions
do ie = nets,nete
do k=1,nlevp; do j=1,np; do i=1,np

do k=1,nlev
dp(:,:, k) =-( (hvcoord%etai(k)-hvcoord%etai(k+1)) * p0)
end do
! no surface topography
p_i(i,j,k) = p0*hvcoord%etai(k)

do k=1,nlevp; do j=1,np; do i=1,np
lon = elem(ie)%spherep(i,j)%lon
lat = elem(ie)%spherep(i,j)%lat


p_i(i,j,k) = p0*hvcoord%etai(k)
call baroclinic_wave_test(is_deep,moist,pertt,dcmip_X,lon,lat,p_i(i,j,k),&
w_i(i,j,k) = 0.0d0
! call this only to compute z_i, will ignore all other output
call baroclinic_wave_test(is_deep,moist,pertt,dcmip_X,lon,lat,p_i(i,j,k),&
z_i(i,j,k),use_zcoords,u(i,j,1),v(i,j,1),T(i,j,1),thetav(i,j,1),phis(i,j),ps(i,j),rho(i,j,1),q(i,j,1,1))
if (k > 1) then
dphi_tmp(i,j,k-1) = (phi_from_z(z_i(i,j,k-1)) - phi_from_z(z_i(i,j,k)))/r_hat_from_phi(phi_from_z(z_i(i,j,k)))**2
end if
end do; end do; end do
q(:,:,:,1:5) = 0.0d0
q(:,:,:,6) = 1
w(:,:,:) = 0.0d0
w_i(:,:,:) = 0.0d0


dphi_tmp = 10

phi_prev(:,:,nlevp) = phis
z_i(:, :, nlevp) = z_from_phi(phis)
pnh_prev(:,:,nlevp) = ps
do k=nlev,1,-1; do j=1,np; do i=1,np

lon = elem(ie)%spherep(i,j)%lon
lat = elem(ie)%spherep(i,j)%lat


! ===================================
if (k.eq. nlev ) then
dmass_i(i,j,k) = -dp(i,j,k)/2
else
dmass_i(i,j,k) = -(dp(i,j,k)+dp(i,j,k+1))/2
end if
dmass(i,j,k) = -dp(i,j,k)

phi_guess(i,j,k) = phi_prev(i,j,k+1) + dphi_tmp(i,j,k)

if (mu_instead_of_eos_rootfinding) then

do newton_idx=1, max_newton
z_guess(i,j,k) = z_from_phi((phi_prev(i,j,k+1) + phi_guess(i,j,k))/2)
z_guess_pert(i,j,k) = z_from_phi((phi_prev(i,j,k+1) +phi_guess(i,j,k)+newton_eps)/2)
call baroclinic_wave_test(is_deep,moist,pertt,dcmip_X,lon,lat, &
pnh_next(i,j,k),z_guess(i,j,k),1,u(i,j,k),v(i,j,k),T(i,j,k),thetav(i,j,k),phis(i,j),ps(i,j),rho(i,j,k),q(i,j,k,1))
call baroclinic_wave_test(is_deep,moist,pertt,dcmip_X,lon,lat, &
pnh_next_pert(i,j,k),z_guess_pert(i,j,k),1,u_pert(i,j,k),v(i,j,k),t_pert(i,j,k),thetav(i,j,k),phis(i,j),ps(i,j),rho(i,j,k),q(i,j,k,1))


r_hat(i,j,k) = r_hat_from_phi(( phi_prev(i,j,k+1) + phi_guess(i,j,k))/2)
r_hat_pert(i,j,k) = r_hat_from_phi(( phi_prev(i,j,k+1) + phi_guess(i,j,k)+newton_eps)/2)

pnh_next(i,j,k) = -Rgas * T(i,j,k) * dmass(i,j,k) / (r_hat(i,j,k)**2 * (phi_guess(i,j,k) - phi_prev(i,j,k+1)))
pnh_next_pert(i,j,k) = -Rgas * t_pert(i,j,k) * dmass(i,j,k) / (r_hat_pert(i,j,k)**2 * (phi_guess(i,j,k) + newton_eps - phi_prev(i,j,k+1)))
!end if
residual(i,j,k) = -u(i,j,k)**2/(z_guess(i,j,k)+rearth) - 2 * omega * cos(lat) * u(i,j,k) + &
g_from_phi(phi_guess(i,j,k)) * (1 - (pnh_next(i,j,k)-pnh_prev(i,j,k+1))/ dmass_i(i,j,k) * r_hat_from_phi(phi_prev(i,j,k+1))**2 )
if (abs(residual(i,j,k)) < 1e-10) then
exit
end if
deriv(i,j,k) = ( (-u_pert(i,j,k)**2/(z_guess_pert(i,j,k)+rearth) - 2 * omega * cos(lat) * u_pert(i,j,k) + &
g_from_phi(phi_guess(i,j,k)) * (1 - (pnh_next_pert(i,j,k)-pnh_prev(i,j,k+1))/ dmass_i(i,j,k) * r_hat_from_phi(phi_prev(i,j,k+1))**2 ))-&
(-u(i,j,k)**2/(z_guess(i,j,k)+rearth) - 2 * omega * cos(lat) * u(i,j,k) + &
g_from_phi(phi_guess(i,j,k)) * (1 - (pnh_next(i,j,k)-pnh_prev(i,j,k+1))/ dmass_i(i,j,k) * r_hat_from_phi(phi_prev(i,j,k+1))**2 )))/newton_eps



if (isnan( residual(i,j,k)/deriv(i,j,k)) .or. phi_guess(i,j,k) < phi_prev(i,j,k+1)) then ! .or. pnh_next(i,j,k) > pnh_prev(i,j,k+1)) then
write(iulog,*)'WARNING: Deep atm geopotential rootfinding encountered zero at lev: ', k
write(iulog, *) "phi next: ", phi_guess(i,j,k), "phi_prev: ", phi_prev(i,j,k+1)
write(iulog, *) "pnh_next: ", pnh_next(i,j,k), "pnh_prev: ", pnh_prev(i,j,k+1)
call abort
end if

phi_guess(i,j,k) = phi_guess(i,j,k) - residual(i,j,k)/deriv(i,j,k)


end do

else
do newton_idx=1, max_newton
z_guess(i,j,k) = z_from_phi((phi_guess(i,j,k) + phi_prev(i,j,k+1))/2)
z_guess_pert(i,j,k) = z_from_phi((phi_prev(i,j,k+1)+phi_guess(i,j,k)+newton_eps)/2)
call baroclinic_wave_test(is_deep,moist,pertt,dcmip_X,lon,lat, &
p_tmp,z_guess(i,j,k),1,u(i,j,k),v(i,j,k),T(i,j,k),thetav(i,j,k),phis(i,j),ps(i,j),rho(i,j,k),q(i,j,k,1))
call baroclinic_wave_test(is_deep,moist,pertt,dcmip_X,lon,lat, &
pnh_next_pert(i,j,k),z_guess_pert(i,j,k),1,u(i,j,k),v(i,j,k),t_pert(i,j,k),thetav(i,j,k),phis(i,j),ps(i,j),rho(i,j,k),q(i,j,k,1))

r_hat(i,j,k) = r_hat_from_phi((phi_guess(i,j,k) + phi_prev(i,j,k+1))/2)
r_hat_pert(i,j,k) = r_hat_from_phi((phi_prev(i,j,k+1) + phi_guess(i,j,k)+newton_eps)/2)
enddo; enddo; enddo
do k=1,nlev; do j=1,np; do i=1,np

! no surface topography
p(i,j,k) = p0*hvcoord%etam(k)
dp(i,j,k) = (hvcoord%etai(k+1)-hvcoord%etai(k))*p0

lon = elem(ie)%spherep(i,j)%lon
lat = elem(ie)%spherep(i,j)%lat

pnh_next(i,j,k) = -Rgas * T(i,j,k) * dmass(i,j,k) / (r_hat(i,j,k)**2 * (phi_guess(i,j,k) - phi_prev(i,j,k+1)))
pnh_next_pert(i,j,k) = -Rgas * t_pert(i,j,k) * dmass(i,j,k) / (r_hat_pert(i,j,k)**2 * (phi_guess(i,j,k) + newton_eps - phi_prev(i,j,k+1)))
residual(i,j,k) = pnh_next(i,j,k) - p_tmp
if (abs(residual(i,j,k)) < 1e-10) then
exit
end if
deriv(i,j,k) = (pnh_next_pert(i,j,k) &
- pnh_next(i,j,k))/newton_eps

if (isnan( residual(i,j,k)/deriv(i,j,k)) .or. phi_guess(i,j,k) < phi_prev(i,j,k+1)) then
write(iulog,*)'WARNING: Deep atm geopotential rootfinding encountered zero at lev: ', k
write(iulog, *) "phi next: ", phi_guess(i,j,k), "phi_prev: ", phi_prev(i,j,k+1)
write(iulog, *) "pnh_next: ", pnh_next(i,j,k), "pnh_prev: ", pnh_prev(i,j,k+1)
call abort
end if

phi_guess(i,j,k) = phi_guess(i,j,k) - residual(i,j,k)/deriv(i,j,k)


end do
q(i,j,k,1:5) = 0.0d0
q(i,j,k,6) = 1
w(i,j,k) = 0.0d0

call baroclinic_wave_test(is_deep,moist,pertt,dcmip_X,lon,lat,p(i,j,k),&
z(i,j,k),use_zcoords,u(i,j,k),v(i,j,k),T(i,j,k),thetav(i,j,k),phis(i,j),ps(i,j),rho(i,j,k),q(i,j,k,1))

! initialize tracer chemistry
call initial_value_terminator( lat*rad2dg, lon*rad2dg, q(i,j,k,4), q(i,j,k,5) )
call set_tracers(q(i,j,k,1:6),6,dp(i,j,k),i,j,k,lat,lon,elem(ie))

end if
if (abs(residual(i,j,k)) > 1e-9) then
write(*,*)'WARNING: Deep atm geopotential rootfinding failed at level ', k, ' residual: ', abs(residual(i,j,k))
call abort
end if
pnh_prev(i,j,k) = pnh_next(i,j,k)
p(i,j,k) = pnh_next(i,j,k)
z_i(i, j, k) = z_from_phi(phi_guess(i, j,k))
z(i,j,k) = z_from_phi((phi_guess(i,j,k) + phi_prev(i,j,k+1))/2)

phi_prev(i,j,k) = phi_guess(i,j,k)
! ===================================
min_thetav = min( min_thetav, thetav(i,j,k) )
max_thetav = max( max_thetav, thetav(i,j,k) )

end do; end do; end do
!do k=1,nlev
! write(*,*) "k: ", k, " max pnh: ", maxval(p(:,:,k)), "min pnh: ", minval(p(:,:,k))
!end do
enddo; enddo; enddo


! do k=1,nlev; do j=1,np; do i=1,np
! ! initialize tracer chemistry
! call initial_value_terminator( lat*rad2dg, lon*rad2dg, q(i,j,k,4), q(i,j,k,5) )
! call set_tracers(q(i,j,k,1:6),6,dp(i,j,k),z(i,j,k),i,j,k,lat,lon,elem(ie))
!
! min_thetav = min( min_thetav, thetav(i,j,k) )
! max_thetav = max( max_thetav, thetav(i,j,k) )
!
! enddo; enddo; enddo

call set_elem_state(u,v,w,w_i,T,ps,phis,p,dp,z,z_i,g,elem(ie),1,nt,ntQ=1)
call tests_finalize(elem(ie),hvcoord)

call pnh_and_exner_from_eos(hvcoord,elem(ie)%state%vtheta_dp(:,:,:,nt),&
elem(ie)%state%dp3d(:,:,:,nt),elem(ie)%state%phinh_i(:,:,:,nt),pnh,exner,dpnh_dp_i)
! call pnh_and_exner_from_eos(hvcoord,T*dp*((p/p0)**(-kappa)),&
! dp,phi_from_z(z_i, nlevp),pnh,exner,dpnh_dp_i)

do k=1,nlevp
if (maxval(abs(dpnh_dp_i(:,:,k)-1)) > 1e-9) then
write(*,*) "dpnh_dp_i error at lev ", k, ": ", maxval(dpnh_dp_i(:,:,k)-1), "minval: ", minval(dpnh_dp_i(:,:,k)-1)
!call abort
end if
end do

enddo
sample_period = 1800.0 ! sec
Expand Down Expand Up @@ -392,7 +269,7 @@ subroutine dcmip2016_test2(elem,hybrid,hvcoord,nets,nete)
q(2)=0
q(3)=0

call set_tracers(q(:),3,dp(i,j,k),z(i,j,k),i,j,k,lat,lon,elem(ie))
call set_tracers(q(:),3,dp(i,j,k),i,j,k,lat,lon,elem(ie))
enddo; enddo; enddo;

call set_elem_state(u,v,w,w_i,T,ps,phis,p,dp,z,z_i,g,elem(ie),1,nt,ntQ=1)
Expand Down Expand Up @@ -498,7 +375,7 @@ subroutine dcmip2016_test3(elem,hybrid,hvcoord,nets,nete)
q (i,j,k,2)= 0 ! no initial clouds
q (i,j,k,3)= 0 ! no initial rain

call set_tracers(q(i,j,k,:),3,dp(i,j,k),z(i,j,k),i,j,k,lat,lon,elem(ie))
call set_tracers(q(i,j,k,:),3,dp(i,j,k),i,j,k,lat,lon,elem(ie))

enddo; enddo
enddo
Expand Down
Loading