From 93623badacf47b1e87b1bb471000ac851de311ae Mon Sep 17 00:00:00 2001 From: Steven van der Linden Date: Sun, 10 Nov 2024 21:16:39 +0100 Subject: [PATCH 1/2] 20241110: added option for pressure gradient (removed hidden dependency that dpl=0 when lcoriol=false) --- src/modboundary.f90 | 12 ++--- src/modforces.f90 | 19 ++++---- src/modglobal.f90 | 15 +++--- src/modstartup.f90 | 116 +++++++++++++++++++++++++++++--------------- 4 files changed, 100 insertions(+), 62 deletions(-) diff --git a/src/modboundary.f90 b/src/modboundary.f90 index 3e63a809..f8a7bd7b 100644 --- a/src/modboundary.f90 +++ b/src/modboundary.f90 @@ -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) diff --git a/src/modforces.f90 b/src/modforces.f90 index e4b0de73..1d453804 100644 --- a/src/modforces.f90 +++ b/src/modforces.f90 @@ -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 @@ -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) diff --git a/src/modglobal.f90 b/src/modglobal.f90 index 188616aa..5f9117e3 100644 --- a/src/modglobal.f90 +++ b/src/modglobal.f90 @@ -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 @@ -367,13 +367,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 @@ -382,7 +380,6 @@ subroutine initglobal ! Variables write(cexpnr,'(i3.3)') iexpnr - ! Create the physical grid variables allocate(dzf(k1), dzfi(k1)) allocate(dzh(k1), dzhi(k1)) diff --git a/src/modstartup.f90 b/src/modstartup.f90 index 22c43648..fcb63758 100644 --- a/src/modstartup.f90 +++ b/src/modstartup.f90 @@ -428,7 +428,7 @@ 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 + use modglobal, only : itot,jtot, ysize,xsize,dtmax,runtime, startfile,lwarmstart,lcoriol,lpressgrad,eps1, imax,jmax use modmpi, only : myid,nprocx,nprocy,mpierr, MPI_FINALIZE use modtimedep, only : ltimedep @@ -474,7 +474,7 @@ subroutine checkinitvalues if (lwarmstart) then if (startfile == '') stop 'no restartfile set' end if - !isurf + !isurf if (myid == 0) then select case (isurf) case(1) @@ -498,6 +498,11 @@ subroutine checkinitvalues 'WARNING: You selected to use time dependent (ltimedep) and heterogeneous surface conditions (lhetero) at the same time' endif + ! 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 @@ -513,7 +518,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,& @@ -524,7 +530,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 @@ -587,10 +592,10 @@ 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) + 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 @@ -960,7 +965,6 @@ subroutine readinitfiles ! 2.1 read and initialise fields !----------------------------------------------------------------- - if(myid==0)then if (ltestbed) then @@ -978,18 +982,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), & @@ -1001,25 +1006,52 @@ 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 @@ -1039,10 +1071,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 @@ -1773,7 +1807,7 @@ 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, & + e12prof, ug, vg, dpdxl, dpdyl, wfls, dqtdxls, dqtdyls, & dqtdtls, dthlrad, kmax) character(*), intent(in) :: filename real(field_r), intent(out) :: height(:) @@ -1784,6 +1818,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(:) @@ -1813,6 +1849,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, & From 640017098909bb9216e18771b3058b16328fafb0 Mon Sep 17 00:00:00 2001 From: Steven van der Linden Date: Fri, 15 Nov 2024 13:46:59 +0100 Subject: [PATCH 2/2] time-dependent profiles of pressure grad and/or Coriolis --- src/modstartup.f90 | 2 ++ src/modtimedep.f90 | 73 +++++++++++++++++++++++++++++++--------------- 2 files changed, 52 insertions(+), 23 deletions(-) diff --git a/src/modstartup.f90 b/src/modstartup.f90 index fcb63758..b41a94d6 100644 --- a/src/modstartup.f90 +++ b/src/modstartup.f90 @@ -1059,6 +1059,8 @@ subroutine readinitfiles 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) diff --git a/src/modtimedep.f90 b/src/modtimedep.f90 index 41948038..befd475f 100644 --- a/src/modtimedep.f90 +++ b/src/modtimedep.f90 @@ -56,6 +56,8 @@ module modtimedep real, allocatable :: timels (:) real, allocatable :: ugt (:,:) real, allocatable :: vgt (:,:) + real, allocatable :: dpdxlt (:,:) + real, allocatable :: dpdylt (:,:) real, allocatable :: wflst (:,:) real, allocatable :: dqtdxlst(:,:) real, allocatable :: dqtdylst(:,:) @@ -67,13 +69,11 @@ module modtimedep real, allocatable :: thlproft(:,:) real, allocatable :: qtproft (:,:) - - contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine inittimedep use modmpi, only :myid,mpierr,comm3d,D_MPI_BCAST - use modglobal, only :cexpnr,k1,kmax,ifinput,runtime,zf,ntimedep + use modglobal, only :cexpnr,k1,kmax,ifinput,runtime,zf,ntimedep,lcoriol,lpressgrad use modsurfdata,only :ps,qts,wqsurf,wtsurf,thls, Qnetav use modtimedepsv, only : inittimedepsv @@ -112,6 +112,8 @@ subroutine inittimedep allocate(timels (0:kls)) allocate(ugt (k1,kls)) allocate(vgt (k1,kls)) + allocate(dpdxlt (k1,kls)) + allocate(dpdylt (k1,kls)) allocate(wflst (k1,kls)) allocate(dqtdxlst (k1,kls)) @@ -139,6 +141,8 @@ subroutine inittimedep ugt = 0 vgt = 0 + dpdxlt = 0 + dpdylt = 0 wflst = 0 dqtdxlst = 0 @@ -245,6 +249,8 @@ subroutine inittimedep stop 'STOP: No time dependend data for end of run' end if end do + + if (ltimedepuv) then ! new, optional format with u,v in ls_flux.inp.* do k=1,kmax @@ -261,26 +267,42 @@ subroutine inittimedep dvdtlst (k,t) end do else - ! old format without u,v in ls_flux.inp.* (default) - do k=1,kmax + ! SvdL, 20241110: if lcoriol, read in 2nd and 3rd column as ug and vg + if (lcoriol) then + ! old format without u,v in ls_flux.inp.* (default) + do k=1,kmax read (ifinput,*) & - height (k) , & - ugt (k,t), & - vgt (k,t), & - wflst (k,t), & - dqtdxlst(k,t), & - dqtdylst(k,t), & - dqtdtlst(k,t), & - thlpcart(k,t) - end do + height (k) , & + ugt (k,t), & + vgt (k,t), & + wflst (k,t), & + dqtdxlst(k,t), & + dqtdylst(k,t), & + dqtdtlst(k,t), & + thlpcart(k,t) + end do + else ! else read in same columns as pressure gradients dpdx and dpdy (chosen for this approach as these columns anyway MUST exist in ls_flux.inp.xxx) + do k=1,kmax + read (ifinput,*) & + height (k) , & + dpdxlt (k,t), & + dpdylt (k,t), & + wflst (k,t), & + dqtdxlst(k,t), & + dqtdylst(k,t), & + dqtdtlst(k,t), & + thlpcart(k,t) + end do + end if + + end if - end do + end do close(ifinput) end if !ltestbed - ! do k=kmax,1,-1 ! write (6,'(3f7.1,5e12.4)') & ! height (k) , & @@ -293,7 +315,6 @@ subroutine inittimedep ! thlpcart(k,t) ! end do - if(timeflux(1)>runtime) then write(6,*) 'Time dependent surface variables do not change before end of' write(6,*) 'simulation. --> only large scale forcings' @@ -320,6 +341,8 @@ subroutine inittimedep call D_MPI_BCAST(timels(1:kls) ,kls ,0,comm3d,mpierr) call D_MPI_BCAST(ugt ,kmax*kls,0,comm3d,mpierr) call D_MPI_BCAST(vgt ,kmax*kls,0,comm3d,mpierr) + call D_MPI_BCAST(dpdxlt ,kmax*kls,0,comm3d,mpierr) ! SvdL, 20241111: broadcast time dependent pressure gradients.. + call D_MPI_BCAST(dpdylt ,kmax*kls,0,comm3d,mpierr) call D_MPI_BCAST(wflst ,kmax*kls,0,comm3d,mpierr) call D_MPI_BCAST(dqtdxlst,kmax*kls ,0,comm3d,mpierr) call D_MPI_BCAST(dqtdylst,kmax*kls ,0,comm3d,mpierr) @@ -381,7 +404,7 @@ subroutine timedepz dvdtls,dvdxls,dvdyls, & dpdxl,dpdyl - use modglobal, only : rtimee,om23_gs,dzf,dzh,k1,kmax,llsadv + use modglobal, only : rtimee,om23_gs,dzf,dzh,k1,kmax,llsadv,lcoriol,lpressgrad use modmpi, only : myid @@ -412,11 +435,15 @@ subroutine timedepz dvdtls = dvdtlst (:,t) + fac * ( dvdtlst (:,t+1) - dvdtlst (:,t) ) thlpcar = thlpcart (:,t) + fac * ( thlpcart (:,t+1) - thlpcart (:,t) ) - - do k=1,kmax - dpdxl(k) = om23_gs*vg(k) - dpdyl(k) = -om23_gs*ug(k) - end do + if (lcoriol) then + do k=1,kmax + dpdxl(k) = om23_gs*vg(k) + dpdyl(k) = -om23_gs*ug(k) + end do + else ! if not lcoriol, update regular pressure gradients + dpdxl = dpdxlt (:,t) + fac * ( dpdxlt (:,t+1) - dpdxlt (:,t) ) + dpdyl = dpdylt (:,t) + fac * ( dpdylt (:,t+1) - dpdylt (:,t) ) + end if whls(1) = 0.0 do k=2,kmax