adv_split_w.F90 14.3 KB
Newer Older
1 2 3 4 5 6
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
! !IROUTINE:  adv_split_w - vertical advection of 3D quantities \label{sec-w-split-adv}
!
! !INTERFACE:
7 8
   subroutine adv_split_w(dt,f,fi,hi,adv,ww,      &
                          splitfac,scheme,tag,az, &
9
                          itersmax,nvd)
10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42
!  Note (KK): Keep in sync with interface in advection_3d.F90
!
! !DESCRIPTION:
!
! Executes an advection step in vertical direction. The 1D advection
! equation
!
! \begin{equation}\label{adv_w_step}
! h^n_{i,j,k} c^n_{i,j,k} =
! h^o_{i,j,k} c^o_{i,j,k}
! - \Delta t
! \left(w_{i,j,k}\tilde c^w_{i,j,k}-w_{i,j,k-1}\tilde c^w_{i,j,k-1}\right),
! \end{equation}
!
! is accompanied by an fractional step for the 1D continuity equation
!
! \begin{equation}\label{adv_w_step_h}
! h^n_{i,j,k}  =
! h^o_{i,j,k}
! - \Delta t
! \left(w_{i,j,k}\tilde -w_{i,j,k-1}\right).
! \end{equation}
!
! Here, $n$ and $o$ denote values before and after this operation,
! respectively, $n$ denote intermediate values when other
! 1D advection steps come after this and $o$ denotes intermediate
! values when other 1D advection steps came before this.
!
! The interfacial fluxes $\tilde c^w_{i,j,k}$ are calculated by means of
! monotone and non-monotone schemes which are described in detail in
! section \ref{sec-u-split-adv} on page \pageref{sec-u-split-adv}.
!
! !USES:
43
   use domain, only: imin,imax,jmin,jmax,kmax,ioff,joff
Knut's avatar
Knut committed
44
   use advection, only: adv_interfacial_reconstruction
Knut's avatar
Knut committed
45
   use advection, only: NOADV,UPSTREAM
46
   use advection_3d, only: W_TAG
47 48 49 50 51
   use halo_zones, only: U_TAG,V_TAG
!$ use omp_lib
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
Knut's avatar
Knut committed
52 53 54 55 56
   REALTYPE,intent(in)                               :: dt,splitfac
   REALTYPE,dimension(I3DFIELD),intent(in),target    :: f
   REALTYPE,dimension(I3DFIELD),intent(in)           :: ww
   integer,intent(in)                                :: scheme,tag,itersmax
   integer,dimension(E2DFIELD),intent(in)            :: az
57 58
!
! !INPUT/OUTPUT PARAMETERS:
59
   REALTYPE,dimension(I3DFIELD),target,intent(inout) :: fi,hi,adv
60
   REALTYPE,dimension(:,:,:),pointer,intent(inout)   :: nvd
61 62
!
! !LOCAL VARIABLES:
63
   logical            :: iterate,use_limiter,calc_nvd,allocated_aux
64
   integer            :: i,j,k,kshift,it,iters,iters_new,rc
65 66
   REALTYPE           :: itersm1,dti,dtik,fio,hio,advn,adv2n,fuu,fu,fd,splitfack
   REALTYPE,dimension(:),allocatable        :: wflux,wflux2
67
   REALTYPE,dimension(:),allocatable,target :: cfl0
68 69
   REALTYPE,dimension(:),pointer            :: fo,faux,fiaux,hiaux,advaux,nvdaux,cfls
   REALTYPE,dimension(:),pointer            :: p_fiaux,p_hiaux,p_advaux,p_nvdaux
Knut's avatar
Knut committed
70
   REALTYPE,dimension(:),pointer            :: p1d
71 72 73 74 75 76 77 78 79 80 81 82
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard & Karsten Bolding
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'adv_split_w() # ',Ncall
#endif
#ifdef SLICE_MODEL
83
   j = jmax/2 ! this MUST NOT be changed!!!
84 85
#endif

86
   if (tag .eq. W_TAG) then
87 88 89 90 91 92
      kshift = 1
   else
      kshift = 0
   end if

   use_limiter = .false.
93
   calc_nvd = associated(nvd)
94 95
   dti = splitfac*dt

96
   iterate = (itersmax .gt. 1)
97 98 99 100 101 102 103
   if (.not. iterate) then
#ifdef NO_BAROTROPIC
      stop 'adv_split_w: do enable iterations with compiler option NO_BAROTROPIC'
#endif
   end if

!$OMP PARALLEL DEFAULT(SHARED)                                         &
104
!$OMP          FIRSTPRIVATE(j,use_limiter)                             &
105 106 107
!$OMP          PRIVATE(rc,wflux,wflux2,cfl0,cfls)                      &
!$OMP          PRIVATE(fo,faux,fiaux,hiaux,advaux,nvdaux)              &
!$OMP          PRIVATE(p_fiaux,p_hiaux,p_advaux,p_nvdaux)              &
108
!$OMP          PRIVATE(itersm1,dtik,splitfack)                         &
109
!$OMP          PRIVATE(i,k,it,iters,iters_new,fio,hio,advn,adv2n,fuu,fu,fd)
110

111 112 113 114 115 116 117

   if (scheme .ne. NOADV) then

!     Each thread allocates its own HEAP storage:
      allocate(wflux(0:kmax),stat=rc)    ! work array
      if (rc /= 0) stop 'adv_split_w: Error allocating memory (wflux)'

118 119 120 121 122
      if (calc_nvd) then
         allocate(wflux2(0:kmax),stat=rc)    ! work array
         if (rc /= 0) stop 'adv_split_w: Error allocating memory (wflux2)'
      end if

123 124 125 126 127 128
#ifndef _POINTER_REMAP_
      allocate(fo(0:kmax),stat=rc)    ! work array
      if (rc /= 0) stop 'adv_split_w: Error allocating memory (fo)'
#endif

      allocated_aux = .false.
129 130 131
#ifdef _POINTER_REMAP_
      if (iterate) then
#endif
132 133 134 135 136
         allocated_aux = .true.

         allocate(fiaux(0:kmax),stat=rc)    ! work array
         if (rc /= 0) stop 'adv_split_w: Error allocating memory (fiaux)'

137 138 139 140 141 142
         allocate(hiaux(0:kmax),stat=rc)    ! work array
         if (rc /= 0) stop 'adv_split_w: Error allocating memory (hiaux)'

         allocate(advaux(0:kmax),stat=rc)    ! work array
         if (rc /= 0) stop 'adv_split_w: Error allocating memory (advaux)'

143 144 145 146 147
         if (calc_nvd) then
            allocate(nvdaux(0:kmax),stat=rc)    ! work array
            if (rc /= 0) stop 'adv_split_w: Error allocating memory (nvdaux)'
         end if

148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163
#ifdef _POINTER_REMAP_
      end if
#endif

      if (scheme.ne.UPSTREAM .or. iterate) then
         allocate(cfl0(0:kmax),stat=rc)    ! work array
         if (rc /= 0) stop 'adv_split_w: Error allocating memory (cfl0)'

         if (iterate) then
            allocate(cfls(0:kmax),stat=rc)    ! work array
            if (rc /= 0) stop 'adv_split_w: Error allocating memory (cfls)'
         else
            cfls => cfl0
         end if
      end if

164
      wflux(0   ) = _ZERO_
165
      wflux(kmax) = _ZERO_
166 167 168 169
      if (calc_nvd) then
         wflux2(0   ) = _ZERO_
         wflux2(kmax) = _ZERO_
      end if
170 171 172 173 174 175 176 177 178 179 180 181

!     Note (KK): as long as h[u|v]n([i|j]max+HALO) are trash (SMALL)
!                they have to be excluded from the loop to avoid
!                unnecessary iterations and warnings

!$OMP DO SCHEDULE(RUNTIME)
#ifndef SLICE_MODEL
      do j=jmin-HALO,jmax+HALO-1
#endif
         do i=imin-HALO,imax+HALO-1
            if (az(i,j) .eq. 1) then
!              Note (KK): exclude vertical advection of normal velocity at open bdys
182 183

               if (scheme.ne.UPSTREAM .or. iterate) then
184
                  do k=1-kshift,kmax-1
185
                     cfl0(k) = abs(ww(i,j,k))*dti/(_HALF_*(hi(i,j,k)+hi(i,j,k+1)))
186
                  end do
187 188 189 190 191 192 193 194
               end if

               iters = 1
               itersm1 = _ONE_
               dtik = dti
               splitfack = splitfac

#ifdef _POINTER_REMAP_
195
               if (.not. iterate) then
Knut's avatar
Knut committed
196 197 198
                  p1d => fi (i,j,:) ; fiaux (0:) => p1d
                  p1d => hi (i,j,:) ; hiaux (0:) => p1d
                  p1d => adv(i,j,:) ; advaux(0:) => p1d
199
                  p1d => nvd(i,j,:) ; nvdaux(0:) => p1d
200 201 202
               end if
#else
               fo = f(i,j,:)
203
#endif
204 205 206 207

               p_fiaux  => fiaux
               p_hiaux  => hiaux
               p_advaux => advaux
208
               p_nvdaux => nvdaux
209 210 211 212 213 214 215 216 217 218 219

               it = 1

               do while (it .le. iters)

                  if (it .eq. 1) then

                     if (allocated_aux) then
                        fiaux  = fi (i,j,:)
                        hiaux  = hi (i,j,:)
                        advaux = adv(i,j,:)
220 221 222
                        if (calc_nvd) then
                           nvdaux = nvd(i,j,:)
                        end if
223
                     end if
224 225

#ifdef _POINTER_REMAP_
Knut's avatar
Knut committed
226
                     p1d => f(i,j,:) ; faux(0:) => p1d
227 228
#else
                     faux => fo
229
#endif
230

231 232 233 234 235
                  end if

                  if (iterate) then
                     if (it .eq. 1) then
                        cfls = cfl0
236
                     else if (scheme.ne.UPSTREAM .or. iters.lt.itersmax) then
237 238 239 240
                        do k=1-kshift,kmax-1
                           cfls(k) = abs(ww(i,j,k))*dti/(_HALF_*(hiaux(k)+hiaux(k+1)))
                        end do
                     end if
241
                     if (iters .lt. itersmax) then
242 243 244
!                       estimate number of iterations by maximum cfl number in water column
                        iters_new = max(1,maxval(ceiling(cfls(1-kshift:kmax-1))))
                        if (iters_new .gt. iters) then
245
                           if (iters_new .gt. itersmax) then
246
!$OMP CRITICAL
247 248
                              STDERR 'adv_split_w: too many iterations needed: ',iters_new
                              STDERR '             at global i j = ',ioff+i,joff+j
249
!$OMP END CRITICAL
250
                              iters = itersmax
251 252 253 254 255 256 257
                           else
                              iters = iters_new
                           end if
                           itersm1 = _ONE_ / iters
                           dtik = dti * itersm1
                           splitfack = splitfac * itersm1
                           if (it .gt. 1) then
258
#ifdef DEBUG
259 260 261 262
!$OMP CRITICAL
                              STDERR 'adv_split_w: restart iterations during it=',it
                              STDERR 'i=',i,' j=',j,':',iters
!$OMP END CRITICAL
263
#endif
264
                              it = 1
265 266 267 268
                              cycle
                           end if
                        end if
                     end if
269
                  end if
270

271 272 273 274
!                 Calculating w-interface fluxes !
                  do k=1-kshift,kmax-1
!                    Note (KK): overwrite zero flux at k=0 in case of W_TAG
                     if (ww(i,j,k) .gt. _ZERO_) then
275
                        fu = faux(k)               ! central
276 277 278 279 280
                        if (scheme .ne. UPSTREAM) then
!                          also fall back to upstream near boundaries
                           use_limiter = (k .gt. 1-kshift)
                        end if
                        if (use_limiter) then
281 282
                           fuu = faux(k-1)            ! upstream
                           fd  = faux(k+1)            ! downstream
283 284
                        end if
                     else
285
                        fu = faux(k+1)               ! central
286 287 288 289 290
                        if (scheme .ne. UPSTREAM) then
!                          also fall back to upstream near boundaries
                           use_limiter = (k .lt. kmax-1)
                        end if
                        if (use_limiter) then
291 292
                           fuu = faux(k+2)            ! upstream
                           fd  = faux(k  )            ! downstream
293 294 295
                        end if
                     end if
                     if (use_limiter) then
Knut's avatar
Knut committed
296
                        fu = adv_interfacial_reconstruction(scheme,cfls(k)*itersm1,fuu,fu,fd)
297
                     end if
298
                     wflux(k) = ww(i,j,k)*fu
299 300 301
                     if (calc_nvd) then
                        wflux2(k) = wflux(k)*fu
                     end if
302
                  end do
303 304 305

#ifdef _POINTER_REMAP_
                  if (it .eq. iters) then
Knut's avatar
Knut committed
306 307 308
                     p1d => fi (i,j,:) ; p_fiaux (0:) => p1d
                     p1d => hi (i,j,:) ; p_hiaux (0:) => p1d
                     p1d => adv(i,j,:) ; p_advaux(0:) => p1d
309 310 311
                     if (calc_nvd) then
                        p1d => nvd(i,j,:) ; p_nvdaux(0:) => p1d
                     end if
312 313 314
                  end if
#endif

315 316
                  do k=1,kmax-kshift
!                    Note (KK): in case of W_TAG do not advect at k=kmax
317
                     fio = fiaux(k)
318
                     hio = hiaux(k)
319
                     p_hiaux(k) = hio - dtik*(ww(i,j,k  )-ww(i,j,k-1))
320
                     advn = splitfack*(wflux(k  )-wflux(k-1))
321
                     p_fiaux(k) = ( hio*fio - dt*advn ) / p_hiaux(k)
322
                     p_advaux(k) = advaux(k) + advn
323 324 325 326 327 328 329 330
                     if (calc_nvd) then
                        adv2n = splitfack*(wflux2(k  )-wflux2(k-1))
!                        p_nvdaux(k) = nvdaux(k)                                            &
!                                     -((p_hiaux(k)*p_fiaux(k)**2 - hio*fio**2)/dt + adv2n)
                        p_nvdaux(k) = ( hio*nvdaux(k)                                        &
                                       -((p_hiaux(k)*p_fiaux(k)**2 - hio*fio**2)/dt + adv2n) &
                                      )/p_hiaux(k)
                     end if
331
                  end do
332 333 334 335

                  faux => p_fiaux
                  it = it + 1

336
               end do
337

338
#ifndef _POINTER_REMAP_
339 340 341
               fi (i,j,1:kmax-kshift) = p_fiaux (1:kmax-kshift)
               hi (i,j,1:kmax-kshift) = p_hiaux (1:kmax-kshift)
               adv(i,j,1:kmax-kshift) = p_advaux(1:kmax-kshift)
342 343 344
               if (calc_nvd) then
                  nvd(i,j,1:kmax-kshift) = p_nvdaux(1:kmax-kshift)
               end if
345
#endif
346

347
            end if
348

349 350 351 352
         end do
#ifndef SLICE_MODEL
      end do
#endif
353
!$OMP END DO NOWAIT
354 355 356 357 358

!     Each thread must deallocate its own HEAP storage:
      deallocate(wflux,stat=rc)
      if (rc /= 0) stop 'adv_split_w: Error deallocating memory (wflux)'

359 360 361 362 363
      if (calc_nvd) then
         deallocate(wflux2,stat=rc)
         if (rc /= 0) stop 'adv_split_w: Error deallocating memory (wflux2)'
      end if

364 365
#ifndef _POINTER_REMAP_
      deallocate(fo,stat=rc)
Knut's avatar
Knut committed
366
      if (rc /= 0) stop 'adv_split_w: Error deallocating memory (fo)'
367 368
#endif

369 370 371
#ifdef _POINTER_REMAP_
      if (iterate) then
#endif
372 373 374 375

         deallocate(fiaux,stat=rc)
         if (rc /= 0) stop 'adv_split_w: Error deallocating memory (fiaux)'

376 377
         deallocate(hiaux,stat=rc)
         if (rc /= 0) stop 'adv_split_w: Error deallocating memory (hiaux)'
378

379 380
         deallocate(advaux,stat=rc)
         if (rc /= 0) stop 'adv_split_w: Error deallocating memory (advaux)'
381

382 383 384 385 386
         if (calc_nvd) then
            deallocate(nvdaux,stat=rc)
            if (rc /= 0) stop 'adv_split_w: Error deallocating memory (nvdaux)'
         end if

387 388 389
#ifdef _POINTER_REMAP_
      end if
#endif
390

391 392 393 394 395 396 397 398 399 400
      if (scheme.ne.UPSTREAM .or. iterate) then
         deallocate(cfl0,stat=rc)
         if (rc /= 0) stop 'adv_split_w: Error deallocating memory (cfl0)'

         if (iterate) then
            deallocate(cfls,stat=rc)
            if (rc /= 0) stop 'adv_split_w: Error deallocating memory (cfls)'
         end if
      end if

401 402 403 404 405 406 407 408 409 410 411 412 413 414
   end if

!$OMP END PARALLEL

#ifdef DEBUG
   write(debug,*) 'Leaving adv_split_w()'
   write(debug,*)
#endif
   return
   end subroutine adv_split_w
!EOC
!-----------------------------------------------------------------------
! Copyright (C) 2004 - Hans Burchard and Karsten Bolding               !
!-----------------------------------------------------------------------