Merge branch 'dev_pressgrad' into dev_pressurefix

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)

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)

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))

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,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)
@@ -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, &