Skip to content
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

Dev pressurefix; horizontal pressure gradient #156

Open
wants to merge 3 commits into
base: dev
Choose a base branch
from
Open

Dev pressurefix; horizontal pressure gradient #156

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
93623ba
20241110: added option for pressure gradient (removed hidden dependen…
Nov 10, 2024
6400170
time-dependent profiles of pressure grad and/or Coriolis
Nov 15, 2024
1d8c33e
Merge branch 'dev_pressgrad' into dev_pressurefix
Feb 7, 2025
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
12 changes: 6 additions & 6 deletions src/modboundary.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -220,12 +220,12 @@ subroutine grwdamp
end do
!$acc end kernels
if(lcoriol) then
!$acc kernels default(present) async(1)
do k=ksp,kmax
up(:,:,k) = up(:,:,k)-(u0(:,:,k)-(ug(k)-cu))*((1./(geodamptime*rnu0))*tsc(k))
vp(:,:,k) = vp(:,:,k)-(v0(:,:,k)-(vg(k)-cv))*((1./(geodamptime*rnu0))*tsc(k))
end do
!$acc end kernels
!$acc kernels default(present) async(1)
do k=ksp,kmax
up(:,:,k) = up(:,:,k)-(u0(:,:,k)-(ug(k)-cu))*((1./(geodamptime*rnu0))*tsc(k))
vp(:,:,k) = vp(:,:,k)-(v0(:,:,k)-(vg(k)-cv))*((1./(geodamptime*rnu0))*tsc(k))
end do
!$acc end kernels
end if
case(2)
!$acc kernels default(present) async(1)
Expand Down
19 changes: 10 additions & 9 deletions src/modforces.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ subroutine forces
! |
!-----------------------------------------------------------------|

use modglobal, only : kmax,dzh,dzf,grav, lpressgrad
use modglobal, only : kmax,dzh,dzf,grav, lpressgrad, lcoriol
use modfields, only : sv0,up,vp,wp,thv0h,dpdxl,dpdyl,thvh
use moduser, only : force_user
use modmicrodata, only : imicro, imicro_bulk, imicro_bin, imicro_sice, imicro_sice2, iqr
Expand All @@ -75,14 +75,15 @@ subroutine forces

if (lforce_user) call force_user

if (lpressgrad) then
!$acc kernels default(present) async(1)
do k = 1, kmax
up(:,:,k) = up(:,:,k) - dpdxl(k) !RN LS pressure gradient force in x,y directions;
vp(:,:,k) = vp(:,:,k) - dpdyl(k)
end do
!$acc end kernels
end if
! SvdL, 20241110: either apply pressure gradient calculated from geostrophic wind speeds (lcoriol) or imposed pressure gradient (lpressgrad)
if (lcoriol .or. lpressgrad) then
!$acc kernels default(present) async(1)
do k = 1, kmax
up(:,:,k) = up(:,:,k) - dpdxl(k) !RN LS pressure gradient force in x,y directions;
vp(:,:,k) = vp(:,:,k) - dpdyl(k)
end do
!$acc end kernels
end if

if((imicro==imicro_sice).or.(imicro==imicro_sice2).or.(imicro==imicro_bulk).or.(imicro==imicro_bin)) then
!$acc kernels default(present) async(2)
Expand Down
15 changes: 6 additions & 9 deletions src/modglobal.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,7 +102,7 @@ module modglobal
real,parameter :: gamma_T_water = 0.57 !< Heat conductivity water [J s-1 m-1 K-1]

logical :: lcoriol = .true. !< switch for coriolis force
logical :: lpressgrad = .true. !< switch for horizontal pressure gradient force
logical :: lpressgrad = .false. !< switch for horizontal pressure gradient force (update 20241110: only intended for channel-like pressure-driven flow, not to be used in combination with coriolis force; thus default to false)
integer :: igrw_damp = 2 !< switch to enable gravity wave damping
real(field_r) :: geodamptime = 7200. !< time scale for nudging to geowind in sponge layer, prevents oscillations
real(field_r) :: uvdamprate = 0. !< rate for damping mean horizontal wind
Expand Down Expand Up @@ -371,13 +371,11 @@ subroutine initglobal
phi = xlat*pi/180.
colat = cos(phi)
silat = sin(phi)
if (lcoriol) then
omega = 7.292e-5
omega_gs = 7.292e-5
else
omega = 0.
omega_gs = 0.
end if

! SvdL, 20241110: define omega(_gs) per default, remove "hidden" rewrite to zero pressure gradient when lcoriol = .false.
omega = 7.292e-5
omega_gs = 7.292e-5

om22 = 2.*omega*colat
om23 = 2.*omega*silat
om22_gs = 2.*omega_gs*colat
Expand All @@ -386,7 +384,6 @@ subroutine initglobal
! Variables
write(cexpnr,'(i3.3)') iexpnr


! Create the physical grid variables
allocate(dzf(k1), dzfi(k1))
allocate(dzh(k1), dzhi(k1))
Expand Down
120 changes: 81 additions & 39 deletions src/modstartup.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -422,7 +422,7 @@ end subroutine startup
subroutine checkinitvalues
use modsurfdata, only: wtsurf, wqsurf, ustin, thls, isurf, ps, lhetero
use modglobal, only: itot, jtot, ysize, xsize, dtmax, runtime, &
startfile, lwarmstart, eps1, imax, jmax, ih, jh
startfile, lwarmstart, lcoriol, lpressgrad, eps1, imax, jmax, ih, jh
use modmpi, only: myid, nprocx, nprocy, mpierr, MPI_FINALIZE
use modtimedep, only: ltimedep

Expand Down Expand Up @@ -513,6 +513,11 @@ subroutine checkinitvalues
end if
end if

! SvdL, 20241110: new lpressgrad intended for pressure-driven channel flow, not to be used i.c.w. Coriolis force
if (lcoriol .and. lpressgrad) then
if (myid==0) stop "Coriolis force (lcoriol) and channel-like pressure gradient (lpressgrad) are mutually exclusive. To use Coriolis force with NO pressure gradient, set geowinds to zero."
end if

end subroutine checkinitvalues

subroutine readinitfiles
Expand All @@ -528,7 +533,8 @@ subroutine readinitfiles
zf,dzf,dzh,rv,rd,cp,rlv,pref0,om23_gs,&
ijtot,cu,cv,e12min,dzh,cexpnr,ifinput,lwarmstart,ltotruntime,itrestart,&
trestart, ladaptive,llsadv,tnextrestart,longint,lconstexner,lopenbc, linithetero, &
lstart_netcdf
lstart_netcdf, &
lcoriol, lpressgrad
use modsubgrid, only : ekm,ekh
use modsurfdata, only : wsvsurf, &
thls,tskin,tskinm,tsoil,tsoilm,phiw,phiwm,Wl,Wlm,thvs,qts,isurf,svs,obl,oblav,&
Expand All @@ -539,7 +545,6 @@ subroutine readinitfiles
use modmpi, only : slabsum,myid,comm3d,mpierr,D_MPI_BCAST
use modthermodynamics, only : thermodynamics,calc_halflev
use moduser, only : initsurf_user

use modtestbed, only : ltestbed,tb_ps,tb_thl,tb_qt,tb_u,tb_v,tb_w,tb_ug,tb_vg,&
tb_dqtdxls,tb_dqtdyls,tb_qtadv,tb_thladv
use modopenboundary, only : openboundary_ghost,openboundary_readboundary,openboundary_initfields
Expand Down Expand Up @@ -602,11 +607,11 @@ subroutine readinitfiles

ps = tb_ps(1)

else if (lstart_netcdf) then
else if (lstart_netcdf) then ! SvdL, 20241110: added reading in of dpdxl and dpdyl, defaults to zero if not present in netcdf file
call init_from_netcdf('init.'//cexpnr//'.nc', height, uprof, vprof, &
thlprof, qtprof, e12prof, ug, vg, wfls, &
dqtdxls, dqtdyls, dqtdtls, thlpcar, kmax)
call tracer_profs_from_netcdf('tracers.'//cexpnr//'.nc', &
thlprof, qtprof, e12prof, ug, vg, dpdxl, dpdyl, wfls, &
dqtdxls, dqtdyls, dqtdtls, thlpcar, kmax)
call tracer_profs_from_netcdf('tracers.'//cexpnr//'.nc', &
tracer_prop, nsv, svprof(1:kmax,:))
else
open (ifinput,file='prof.inp.'//cexpnr,status='old',iostat=ierr)
Expand Down Expand Up @@ -974,7 +979,6 @@ subroutine readinitfiles
! 2.1 read and initialise fields
!-----------------------------------------------------------------


if(myid==0)then

if (ltestbed) then
Expand All @@ -992,18 +996,19 @@ subroutine readinitfiles
thlpcar(k) = tb_thladv(1,k)
end do

else if (lstart_netcdf) then
continue! Profiles have been read by init_from_netcdf
else
open (ifinput,file='lscale.inp.'//cexpnr, status='old',iostat=ierr)
if (ierr /= 0) then
write(6,*) 'Cannot open the file ', 'lscale.inp.'//cexpnr
STOP
end if
read (ifinput,'(a80)') chmess
read (ifinput,'(a80)') chmess

if (lstart_netcdf) then
continue ! Profiles have been read by init_from_netcdf
else
open (ifinput,file='lscale.inp.'//cexpnr, status='old',iostat=ierr)
if (ierr /= 0) then
write(6,*) 'Cannot open the file ', 'lscale.inp.'//cexpnr
STOP
end if
read (ifinput,'(a80)') chmess
read (ifinput,'(a80)') chmess
! SvdL, 20241110: if lcoriol, read in 2nd and 3rd column as ug and vg
if (lcoriol) then
do k=1,kmax
read (ifinput,*) &
height (k), &
Expand All @@ -1015,32 +1020,61 @@ subroutine readinitfiles
dqtdtls(k), &
thlpcar(k)
end do
close(ifinput)
else ! else read in same columns as pressure gradients dpdx and dpdy (chosen for this approach as these columns anyway MUST exist in lscale.inp.xxx)
do k=1,kmax
read (ifinput,*) &
height (k), &
dpdxl (k), &
dpdyl (k), &
wfls (k), &
dqtdxls(k), &
dqtdyls(k), &
dqtdtls(k), &
thlpcar(k)
end do
end if

close(ifinput)
end if

write(6,*) ' height u_geo v_geo subs ' &
,' dqtdx dqtdy dqtdtls thl_rad '
do k=kmax,1,-1
write (6,'(3f7.1,5e12.4)') &
zf (k), &
ug (k), &
vg (k), &
wfls (k), &
dqtdxls(k), &
dqtdyls(k), &
dqtdtls(k), &
thlpcar(k)
end do

if (.not. lpressgrad) then
write(6,*) ' height u_geo v_geo subs ' &
,' dqtdx dqtdy dqtdtls thl_rad '
do k=kmax,1,-1
write (6,'(3f7.1,5e12.4)') &
zf (k), &
ug (k), &
vg (k), &
wfls (k), &
dqtdxls(k), &
dqtdyls(k), &
dqtdtls(k), &
thlpcar(k)
end do
else
write(6,*) ' height dpdxl dpdyl subs ' &
,' dqtdx dqtdy dqtdtls thl_rad '
do k=kmax,1,-1
write (6,'(3f7.1,5e12.4)') &
zf (k), &
dpdxl (k), &
dpdyl (k), &
wfls (k), &
dqtdxls(k), &
dqtdyls(k), &
dqtdtls(k), &
thlpcar(k)
end do
end if

end if ! end myid==0

! MPI broadcast variables read in

call D_MPI_BCAST(ug ,kmax,0,comm3d,mpierr)
call D_MPI_BCAST(vg ,kmax,0,comm3d,mpierr)
call D_MPI_BCAST(dpdxl ,kmax,0,comm3d,mpierr) ! SvdL, 20241111: broadcast new pressure gradients
call D_MPI_BCAST(dpdyl ,kmax,0,comm3d,mpierr)
call D_MPI_BCAST(wfls ,kmax,0,comm3d,mpierr)
call D_MPI_BCAST(dqtdxls ,kmax,0,comm3d,mpierr)
call D_MPI_BCAST(dqtdyls ,kmax,0,comm3d,mpierr)
Expand All @@ -1053,10 +1087,12 @@ subroutine readinitfiles

!******include rho if rho = rho(z) /= 1.0 ***********

do k = 1, kmax
dpdxl(k) = om23_gs*vg(k)
dpdyl(k) = -om23_gs*ug(k)
end do
if (lcoriol) then ! only when Coriolis is enabled, calculte pressure gradients via geostrophic wind speeds (otherwise assume they have been set already)
do k = 1, kmax
dpdxl(k) = om23_gs*vg(k)
dpdyl(k) = -om23_gs*ug(k)
end do
end if

!-----------------------------------------------------------------
! 2.5 make large-scale horizontal gradients
Expand Down Expand Up @@ -1787,8 +1823,8 @@ end subroutine baseprofs
!! \note Tracers are read from tracers.XXX.nc, not here.
!! \todo Make DEPHY-compatible.
subroutine init_from_netcdf(filename, height, uprof, vprof, thlprof, qtprof, &
e12prof, ug, vg, wfls, dqtdxls, dqtdyls, &
dqtdtls, dthlrad, kmax)
e12prof, ug, vg, dpdxl, dpdyl, wfls, dqtdxls, dqtdyls, &
dqtdtls, dthlrad, kmax)
character(*), intent(in) :: filename
real(field_r), intent(out) :: height(:)
real(field_r), intent(out) :: uprof(:)
Expand All @@ -1798,6 +1834,8 @@ subroutine init_from_netcdf(filename, height, uprof, vprof, thlprof, qtprof, &
real(field_r), intent(out) :: e12prof(:)
real(field_r), intent(out) :: ug(:)
real(field_r), intent(out) :: vg(:)
real(field_r), intent(out) :: dpdxl(:)
real(field_r), intent(out) :: dpdyl(:)
real(field_r), intent(out) :: wfls(:)
real(field_r), intent(out) :: dqtdxls(:)
real(field_r), intent(out) :: dqtdyls(:)
Expand Down Expand Up @@ -1827,6 +1865,10 @@ subroutine init_from_netcdf(filename, height, uprof, vprof, thlprof, qtprof, &
fillvalue=0._field_r)
call read_nc_field(ncid, "vg", vg, start=1, count=kmax, &
fillvalue=0._field_r)
call read_nc_field(ncid, "dpdx", dpdxl, start=1, count=kmax, &
fillvalue=0._field_r)
call read_nc_field(ncid, "dpdy", dpdyl, start=1, count=kmax, &
fillvalue=0._field_r)
call read_nc_field(ncid, "wa", wfls, start=1, count=kmax, &
fillvalue=0._field_r)
call read_nc_field(ncid, "dqtdxls", dqtdxls, start=1, count=kmax, &
Expand Down
Loading