Merge pull request #4 from CasparJungbacker/OpenACC-subgrid

Offload subgrid

.github/workflows/precision.yml

@@ -4,7 +4,13 @@ name: DALES precision CI
 on:
   # Triggers the workflow on push or pull request events but only for the feature/single_precision branch
   push:
+    branches:
+      - master
+      - OpenACC-Dev
   pull_request:
+    branches:
+      - master
+      - OpenACC-Dev
   # Allows you to run this workflow manually from the Actions tab
   workflow_dispatch:
 

src/modadvection.f90

@@ -47,7 +47,7 @@ subroutine advection
 
   ! leq = .false. ! for testing that the non-uniform advection routines agree with the uniform ones
                   ! when the grid is uniform
-  !$acc data copyin(u0, v0, w0, e120, thl0, qt0, dzf, dzh, rhobf, rhobh) copy(up, vp, wp, e12p, thlp, qtp)
+  !$acc enter data copyin(u0, v0, w0, e120, thl0, qt0, dzf, dzh, rhobf, rhobh, up, vp, wp, e12p, thlp, qtp)
   select case(iadv_mom)
     case(iadv_cd2)
       call advecu_2nd(u0,up)
@@ -212,6 +212,5 @@ subroutine advection
     end select
   end do
   !$acc wait
-  !$acc end data
 end subroutine advection
 end module modadvection

src/modbudget.f90

@@ -213,7 +213,7 @@ subroutine do_genbudget
                           ijtot,cu,cv,grav, &
                           dxi,dyi,dx2i,dy2i
     use modsurfdata,only : ustar
-    use modsubgriddata, only : ekm
+    use modsubgriddata, only : ekm,anis_fac
     use modpois,    only : p
     use modfields,  only : u0,v0,w0,thv0h,u0av,v0av,rhobf,rhobh,thvh
 !cstep    use modtilt,    only : adjustbudget,ltilted
@@ -407,12 +407,12 @@ subroutine do_genbudget
 
       tau1j(i,j,k) = &
                ( ekm(i,j,k)  * (u0(i+1,j,k)-u0(i,j,k)) &
-                -ekm(i-1,j,k)* (u0(i,j,k)-u0(i-1,j,k)) ) * 2. * dx2i &
+                -ekm(i-1,j,k)* (u0(i,j,k)-u0(i-1,j,k)) ) * 2. * dx2i * anis_fac(k) &
                + &
                ( empo * ( (u0(i,jp,k)-u0(i,j,k))   *dyi &
                          +(v0(i,jp,k)-v0(i-1,jp,k))*dxi) &
                 -emmo * ( (u0(i,j,k)-u0(i,jm,k))   *dyi &
-                         +(v0(i,j,k)-v0(i-1,j,k))  *dxi) ) * dyi &
+			+(v0(i,j,k)-v0(i-1,j,k))  *dxi) ) * dyi * anis_fac(k) &
                + &
                ( rhobh(kp)/rhobf(k) * emop * ( (u0(i,j,kp)-u0(i,j,k))   /dzh(kp) &
                          +(w0(i,j,kp)-w0(i-1,j,kp))*dxi) &
@@ -455,10 +455,10 @@ subroutine do_genbudget
                ( epmo * ( (v0(i+1,j,k)-v0(i,j,k))   *dxi &
                          +(u0(i+1,j,k)-u0(i+1,jm,k))*dyi) &
                 -emmo * ( (v0(i,j,k)-v0(i-1,j,k))   *dxi &
-                         +(u0(i,j,k)-u0(i,jm,k))    *dyi)  ) * dxi &
+                         +(u0(i,j,k)-u0(i,jm,k))    *dyi)  ) * dxi * anis_fac(k) &
                + &
                ( ekm(i,j,k) * (v0(i,jp,k)-v0(i,j,k)) &
-                -ekm(i,jm,k)* (v0(i,j,k)-v0(i,jm,k))  ) * 2. * dy2i &
+                -ekm(i,jm,k)* (v0(i,j,k)-v0(i,jm,k))  ) * 2. * dy2i * anis_fac(k) &
                + &
                ( rhobh(kp)/rhobf(k) * eomp * ( (v0(i,j,kp)-v0(i,j,k))    /dzh(kp) &
                          +(w0(i,j,kp)-w0(i,jm,kp))  *dyi) &
@@ -501,12 +501,12 @@ subroutine do_genbudget
                  ( epom * ( (w0(i+1,j,k)-w0(i,j,k))    *dxi &
                            +(u0(i+1,j,k)-u0(i+1,j,km)) /dzh(k) ) &
                   -emom * ( (w0(i,j,k)-w0(i-1,j,k))    *dxi &
-                           +(u0(i,j,k)-u0(i,j,km))     /dzh(k) ))*dxi &
+                           +(u0(i,j,k)-u0(i,j,km))     /dzh(k) ))*dxi * anis_fac(k) &
                  + &
                  ( eopm * ( (w0(i,jp,k)-w0(i,j,k))     *dyi &
                            +(v0(i,jp,k)-v0(i,jp,km))   /dzh(k) ) &
                   -eomm * ( (w0(i,j,k)-w0(i,jm,k))     *dyi &
-                           +(v0(i,j,k)-v0(i,j,km))     /dzh(k) ))*dyi &
+                           +(v0(i,j,k)-v0(i,j,km))     /dzh(k) ))*dyi * anis_fac(k) &
                  + &
                   ( rhobf(k)/rhobh(k) * ekm(i,j,k) * (w0(i,j,kp)-w0(i,j,k)) /dzf(k) &
                   - rhobf(km)/rhobh(k) * ekm(i,j,km)* (w0(i,j,k)-w0(i,j,km)) /dzf(km) ) * 2. &
@@ -553,12 +553,12 @@ subroutine do_genbudget
 
     tau1j(i,j,1) = &
              ( ekm(i,j,1)  * (u0(i+1,j,1)-u0(i,j,1)) &
-              -ekm(i-1,j,1)* (u0(i,j,1)-u0(i-1,j,1)) ) * 2. * dx2i &
+              -ekm(i-1,j,1)* (u0(i,j,1)-u0(i-1,j,1)) ) * 2. * dx2i * anis_fac(k) &
              + &
              ( empo * ( (u0(i,jp,1)-u0(i,j,1))   *dyi &
                        +(v0(i,jp,1)-v0(i-1,jp,1))*dxi) &
              - emmo * ( (u0(i,j,1)-u0(i,jm,1))   *dyi &
-                       +(v0(i,j,1)-v0(i-1,j,1))  *dxi)   ) * dyi &
+                       +(v0(i,j,1)-v0(i-1,j,1))  *dxi)   ) * dyi * anis_fac(k) &
              + &
              ( rhobh(2)/rhobf(1) * emop * ( (u0(i,j,2)-u0(i,j,1))    /dzh(2) &
                        +(w0(i,j,2)-w0(i-1,j,2))  *dxi) &
@@ -584,10 +584,10 @@ subroutine do_genbudget
             ( epmo * ( (v0(i+1,j,1)-v0(i,j,1))   *dxi &
                       +(u0(i+1,j,1)-u0(i+1,jm,1))*dyi) &
              -emmo * ( (v0(i,j,1)-v0(i-1,j,1))   *dxi &
-                      +(u0(i,j,1)-u0(i,jm,1))    *dyi)   ) * dxi &
+                      +(u0(i,j,1)-u0(i,jm,1))    *dyi)   ) * dxi * anis_fac(k) &
             + &
             ( ekm(i,j,1) * (v0(i,jp,1)-v0(i,j,1)) &
-             -ekm(i,jm,1)* (v0(i,j,1)-v0(i,jm,1))  ) * 2. * dy2i &
+             -ekm(i,jm,1)* (v0(i,j,1)-v0(i,jm,1))  ) * 2. * dy2i * anis_fac(k) &
             + &
             ( rhobh(2)/rhobf(1) * eomp * ( (v0(i,j,2)-v0(i,j,1))     /dzh(2) &
                       +(w0(i,j,2)-w0(i,jm,2))    *dyi) &

src/modforces.f90

@@ -63,69 +63,36 @@ subroutine forces
   use modglobal, only : i1,j1,kmax,dzh,dzf,grav, lpressgrad
   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,iqr
+  use modmicrodata, only : imicro, imicro_bulk, imicro_bin, imicro_sice, imicro_sice2, iqr
   implicit none
 
-  integer i, j, k, jm, jp, km, kp
+  integer i, j, k
 
   if (lforce_user) call force_user
 
-  if((imicro==imicro_sice).or.(imicro==imicro_bulk).or.(imicro==imicro_bin)) then
+  if (lpressgrad) then
+     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
+  end if
+
+  if((imicro==imicro_sice).or.(imicro==imicro_sice2).or.(imicro==imicro_bulk).or.(imicro==imicro_bin)) then
     do k=2,kmax
-      kp=k+1
-      km=k-1
-    do j=2,j1
-      jp=j+1
-      jm=j-1
-    do i=2,i1
-
-    if (lpressgrad) then
-      up(i,j,k) = up(i,j,k) - dpdxl(k)      !RN LS pressure gradient force in x,y directions;
-      vp(i,j,k) = vp(i,j,k) - dpdyl(k)
-    end if
-      wp(i,j,k) = wp(i,j,k) + grav*(thv0h(i,j,k)-thvh(k))/thvh(k) - &
-                  grav*(sv0(i,j,k,iqr)*dzf(k-1)+sv0(i,j,k-1,iqr)*dzf(k))/(2.0*dzh(k))
-    end do
-    end do
+       wp(:,:,k) = wp(:,:,k) + grav*(thv0h(:,:,k)-thvh(k))/thvh(k) - &
+                  grav*(sv0(:,:,k,iqr)*dzf(k-1)+sv0(:,:,k-1,iqr)*dzf(k))/(2.0*dzh(k))
     end do
   else
     do k=2,kmax
-      kp=k+1
-      km=k-1
-    do j=2,j1
-      jp=j+1
-      jm=j-1
-    do i=2,i1
-      up(i,j,k) = up(i,j,k) - dpdxl(k)
-      vp(i,j,k) = vp(i,j,k) - dpdyl(k)
-      wp(i,j,k) = wp(i,j,k) + grav*(thv0h(i,j,k)-thvh(k))/thvh(k)
-    end do
-    end do
+      wp(:,:,k) = wp(:,:,k) + grav*(thv0h(:,:,k)-thvh(k))/thvh(k)
     end do
   end if
 
 !     --------------------------------------------
 !     special treatment for lowest full level: k=1
 !     --------------------------------------------
+  wp(:,:,1) = 0.0
 
-  do j=2,j1
-    jp = j+1
-    jm = j-1
-  do i=2,i1
-
-    if (lpressgrad) then
-      up(i,j,1) = up(i,j,1) - dpdxl(1)
-      vp(i,j,1) = vp(i,j,1) - dpdyl(1)
-    end if
-
-    wp(i,j,1) = 0.0
-
-  end do
-  end do
-!     ----------------------------------------------end i,j-loop
-
-
-  return
   end subroutine forces
   subroutine coriolis
 
@@ -152,7 +119,7 @@ subroutine coriolis
   integer i, j, k, jm, jp, km, kp
 
   if (lcoriol .eqv. .false.) return
-  
+
   do k=2,kmax
     kp=k+1
     km=k-1
@@ -232,80 +199,45 @@ subroutine lstend
                         dqtdtls, dthldtls, dudtls, dvdtls
   implicit none
 
-  integer i,j,k,n,kp,km
-  real subs_thl,subs_qt,subs_sv,subs_u,subs_v
+  integer k,kp,km
 
 !     1. DETERMINE LARGE SCALE TENDENCIES
 !        --------------------------------
 
 !     1.1 lowest model level above surface : only downward component
-
-  subs_thl = 0.
-  subs_qt  = 0.
-  subs_sv  = 0.
-  subs_u   = 0.
-  subs_v   = 0.
-
-  do j=2,j1
-    do i=2,i1
-      k = 1
-      if (whls(2).lt.0) then !neglect effect of mean ascending on tendencies at the lowest full level
-        subs_thl     = 0.5*whls(2)  *(thl0(i,j,2)-thl0(i,j,1))/dzh(2)
-        subs_qt      = 0.5*whls(2)  *(qt0(i,j,2)-qt0(i,j,1) )/dzh(2)
-        if (lmomsubs) then
-          subs_u     = 0.5*whls(2)  *(u0(i,j,2)-u0(i,j,1))/dzh(2)
-          subs_v     = 0.5*whls(2)  *(v0(i,j,2)-v0(i,j,1))/dzh(2)
-        endif
-        do n=1,nsv
-          subs_sv =  0.5*whls(2)  *(sv0(i,j,2,n)-sv0(i,j,1,n)  )/dzh(2)
-          svp(i,j,1,n) = svp(i,j,1,n)-subs_sv
-        enddo
-      endif
-      thlp(i,j,1) = thlp(i,j,1) -u0av(1)*dthldxls(1)-v0av(1)*dthldyls(1)-subs_thl + dthldtls(1)
-      qtp(i,j,1)  = qtp (i,j,1) -u0av(1)*dqtdxls (1)-v0av(1)*dqtdyls (1)-subs_qt  + dqtdtls(1)
-      up  (i,j,1) = up  (i,j,1) -u0av(1)*dudxls  (1)-v0av(1)*dudyls  (1)-subs_u   + dudtls(1)
-      vp  (i,j,1) = vp  (i,j,1) -u0av(1)*dvdxls  (1)-v0av(1)*dvdyls  (1)-subs_v   + dvdtls(1)
-    end do
-  end do
-
 !     1.2 other model levels twostream
 
-  do k=2,kmax
+  do k=1,kmax
     kp=k+1
     km=k-1
-    do j=2,j1
-      do i=2,i1
-        if (whls(kp).lt.0) then   !downwind scheme for subsidence
-          subs_thl    = whls(kp) * (thl0(i,j,kp) - thl0(i,j,k))/dzh(kp)
-          subs_qt     = whls(kp) * (qt0 (i,j,kp) - qt0 (i,j,k))/dzh(kp)
-          if (lmomsubs) then
-            subs_u    = whls(kp) * (u0(i,j,kp) - u0(i,j,k))/dzh(kp)
-            subs_v    = whls(kp) * (v0(i,j,kp) - v0(i,j,k))/dzh(kp)
-          endif
-          do n=1,nsv
-            subs_sv   = whls(kp)  *(sv0(i,j,kp,n) - sv0(i,j,k,n))/dzh(kp)
-            svp(i,j,k,n) = svp(i,j,k,n)-subs_sv
-          enddo
-        else !downwind scheme for mean upward motions
-          subs_thl    = whls(k) * (thl0(i,j,k) - thl0(i,j,km))/dzh(k)
-          subs_qt     = whls(k) * (qt0 (i,j,k) - qt0 (i,j,km))/dzh(k)
+
+    if (whls(kp).lt.0) then   !downwind scheme for subsidence
+       thlp(2:i1,2:j1,k) = thlp(2:i1,2:j1,k) - whls(kp) * (thl0(2:i1,2:j1,kp) - thl0(2:i1,2:j1,k))/dzh(kp)
+       qtp (2:i1,2:j1,k) = qtp (2:i1,2:j1,k) - whls(kp) * (qt0 (2:i1,2:j1,kp) - qt0 (2:i1,2:j1,k))/dzh(kp)
+       if (lmomsubs) then
+          up(2:i1,2:j1,k) = up(2:i1,2:j1,k) - whls(kp) * (u0(2:i1,2:j1,kp) - u0(2:i1,2:j1,k))/dzh(kp)
+          vp(2:i1,2:j1,k) = vp(2:i1,2:j1,k) - whls(kp) * (v0(2:i1,2:j1,kp) - v0(2:i1,2:j1,k))/dzh(kp)
+       endif
+
+       svp(2:i1,2:j1,k,:) = svp(2:i1,2:j1,k,:) - whls(kp) * (sv0(2:i1,2:j1,kp,:) - sv0(2:i1,2:j1,k,:))/dzh(kp)
+
+    else !downwind scheme for mean upward motions
+       if (k > 1) then !neglect effect of mean ascending on tendencies at the lowest full level
+          thlp(2:i1,2:j1,k) = thlp(2:i1,2:j1,k) - whls(k) * (thl0(2:i1,2:j1,k) - thl0(2:i1,2:j1,km))/dzh(k)
+          qtp (2:i1,2:j1,k) = qtp (2:i1,2:j1,k) - whls(k) * (qt0 (2:i1,2:j1,k) - qt0 (2:i1,2:j1,km))/dzh(k)
           if (lmomsubs) then
-            subs_u    = whls(k) * (u0(i,j,k) - u0(i,j,km))/dzh(k)
-            subs_v    = whls(k) * (v0(i,j,k) - v0(i,j,km))/dzh(k)
+             up(2:i1,2:j1,k) = up(2:i1,2:j1,k) - whls(k) * (u0(2:i1,2:j1,k) - u0(2:i1,2:j1,km))/dzh(k)
+             vp(2:i1,2:j1,k) = vp(2:i1,2:j1,k) - whls(k) * (v0(2:i1,2:j1,k) - v0(2:i1,2:j1,km))/dzh(k)
           endif
-          do n=1,nsv
-            subs_sv   = whls(k) * (sv0(i,j,k,n) - sv0(i,j,km,n))/dzh(k)
-            svp(i,j,k,n) = svp(i,j,k,n)-subs_sv
-          enddo
-        endif
-
-        thlp(i,j,k) = thlp(i,j,k)-u0av(k)*dthldxls(k)-v0av(k)*dthldyls(k)-subs_thl + dthldtls(k)
-        qtp (i,j,k) = qtp (i,j,k)-u0av(k)*dqtdxls (k)-v0av(k)*dqtdyls (k)-subs_qt  + dqtdtls(k)
-        up  (i,j,k) = up  (i,j,k)-u0av(k)*dudxls  (k)-v0av(k)*dudyls  (k)-subs_u   + dudtls(k)
-        vp  (i,j,k) = vp  (i,j,k)-u0av(k)*dvdxls  (k)-v0av(k)*dvdyls  (k)-subs_v   + dvdtls(k)
-
-      enddo
-    enddo
+          svp(2:i1,2:j1,k,:) = svp(2:i1,2:j1,k,:)-whls(k) * (sv0(2:i1,2:j1,k,:) - sv0(2:i1,2:j1,km,:))/dzh(k)
+       endif
+    endif
+
+    thlp(2:i1,2:j1,k) = thlp(2:i1,2:j1,k)-u0av(k)*dthldxls(k)-v0av(k)*dthldyls(k) + dthldtls(k)
+    qtp (2:i1,2:j1,k) = qtp (2:i1,2:j1,k)-u0av(k)*dqtdxls (k)-v0av(k)*dqtdyls (k) + dqtdtls(k)
+    up  (2:i1,2:j1,k) = up  (2:i1,2:j1,k)-u0av(k)*dudxls  (k)-v0av(k)*dudyls  (k) + dudtls(k)
+    vp  (2:i1,2:j1,k) = vp  (2:i1,2:j1,k)-u0av(k)*dvdxls  (k)-v0av(k)*dvdyls  (k) + dvdtls(k)
+
   enddo
 
   return

src/modglobal.f90

@@ -216,7 +216,7 @@ module modglobal
       real :: xsize    = -1 !<  domain size in x-direction
       real :: ysize    = -1 !<  domain size in y-direction
       real, allocatable :: delta(:)       !<  (dx*dy*dz)**(1/3)
-      real, allocatable :: deltai(:)       !<  (dx*dy*dz)**(-1/3)
+      real, allocatable :: deltai(:)       !<  (dx*dy*dz)**(-1/3)  or dzf**-1 for anisotropic diffusion
 
       logical :: leq      = .true.  !<  switch for (non)-equidistant mode.
       logical :: lmomsubs = .false.  !<  switch to apply subsidence on the momentum or not
@@ -433,7 +433,7 @@ subroutine initglobal
     do k=1,k1
 
        delta(k) = (dx*dy*dzf(k))**(1./3.)
-       deltai(k) = 1./delta(k)
+       deltai(k) = 1./delta(k)     !can be overruled in modsubgrid in case anisotropic diffusion is applied
     end do
 
   !--------------------------------------------------