diff --git a/src/modboundary.f90 b/src/modboundary.f90 index 3e63a8096..f8a7bd7be 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 e718aa41e..83f72a636 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 95080af77..7c9e2c82d 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 @@ -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 @@ -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)) diff --git a/src/modstartup.f90 b/src/modstartup.f90 index 44f601a4d..742e47903 100644 --- a/src/modstartup.f90 +++ b/src/modstartup.f90 @@ -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 @@ -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 @@ -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,& @@ -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 @@ -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) @@ -974,7 +979,6 @@ subroutine readinitfiles ! 2.1 read and initialise fields !----------------------------------------------------------------- - if(myid==0)then if (ltestbed) then @@ -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), & @@ -1015,25 +1020,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 @@ -1041,6 +1073,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) @@ -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 @@ -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(:) @@ -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(:) @@ -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, & diff --git a/src/modtimedep.f90 b/src/modtimedep.f90 index 419480381..befd475fb 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