!> \file tstep.f90 !! Performs the time integration !> !! Performs the time integration !> !! Tstep uses adaptive timestepping and 3rd order Runge Kutta time integration. !! The adaptive timestepping chooses it's delta_t according to the courant number !! and the diffusion number, depending on the advection scheme in use. !! !! \author Jasper Tomas, TU Delft !! \author Chiel van Heerwaarden, Wageningen University !! \author Thijs Heus,MPI-M !! \see Wicker and Skamarock 2002 !! \par Revision list ! This file is part of DALES. ! ! DALES is free software; you can redistribute it and/or modify ! it under the terms of the GNU General Public License as published by ! the Free Software Foundation; either version 3 of the License, or ! (at your option) any later version. ! ! DALES is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU General Public License for more details. ! ! You should have received a copy of the GNU General Public License ! along with this program. If not, see <http://www.gnu.org/licenses/>. ! ! Copyright 1993-2009 Delft University of Technology, Wageningen University, Utrecht University, KNMI ! !> Determine time step size dt in initialization and update time variables !! !! The size of the timestep Delta t is determined adaptively, and is limited by both the Courant-Friedrichs-Lewy criterion CFL !! \latexonly !! \begin{equation} !! \CFL = \mr{max}\left(\left|\frac{u_i \Delta t}{\Delta x_i}\right|\right), !! \end{equation} !! and the diffusion number $d$. The timestep is further limited by the needs of other modules, e.g. the statistics. !! \endlatexonly module modtstep implicit none contains subroutine tstep_update use modglobal, only : ib,ie,jb,je,rk3step,timee,dtmax,dt,ntimee,ntrun,courant,diffnr,& kb,ke,dxi,dx2i,dyi,dy2i,dzh,dt_lim,ladaptive,timeleft,dt,lwarmstart,& dzh2i use modfields, only : um,vm,wm use modsubgriddata, only : ekm,ekh use modmpi, only : comm3d,mpierr,my_real,mpi_max implicit none integer :: i, j, k real,save :: courtot=-1.,diffnrtot=-1. real :: courtotl,courold,diffnrtotl,diffnrold ! logical,save :: spinup=.true. logical,save :: spinup=.false. if(lwarmstart) spinup = .false. rk3step = mod(rk3step,3) + 1 if(rk3step == 1) then ! Initialization if (spinup) then write(6,*) '!spinup!' if (ladaptive) then courold = courtot diffnrold = diffnrtot courtotl=0. diffnrtotl = 0. do k=kb,ke do j=jb,je do i=ib,ie courtotl = max(courtotl,(abs(um(i,j,k))*dxi + abs(vm(i,j,k))*dyi + abs(wm(i,j,k))/dzh(k))*dt) ! diffnrtotl = max(diffnrtotl, ekm(i,j,k)*(1/dzh(k)**2 + dxh2i(i) + dy2i)*dt ) diffnrtotl = max(diffnrtotl, ekm(i,j,k)*(dzh2i(k) + dx2i + dy2i)*dt, & ekh(i,j,k)*(dzh2i(k) + dx2i + dy2i)*dt ) end do end do end do call MPI_ALLREDUCE(courtotl,courtot,1,MY_REAL,MPI_MAX,comm3d,mpierr) call MPI_ALLREDUCE(diffnrtotl,diffnrtot,1,MY_REAL,MPI_MAX,comm3d,mpierr) if ( diffnrold>0) then dt = min(dtmax,dt*courant/courtot,dt*diffnr/diffnrtot) if ((abs(courtot-courold)/courold<0.1) .and. (abs(diffnrtot-diffnrold)/diffnrold<0.1)) then spinup = .false. end if end if dt = dt dt_lim = timeleft timee = timee + dt timeleft = timeleft- dt ntimee = ntimee + 1 ntrun = ntrun + 1 else dt = 2 * dt if (dt >= dtmax) then dt = dtmax spinup = .false. end if end if ! Normal time loop else !spinup = .false. if (ladaptive) then courtotl=0. diffnrtotl = 1e-5 do k=kb,ke do j=jb,je do i=ib,ie courtotl = max(courtotl,(abs(um(i,j,k))*dxi + abs(vm(i,j,k))*dyi + abs(wm(i,j,k))/dzh(k))*dt) diffnrtotl = max(diffnrtotl, ekm(i,j,k)*(dzh2i(k) + dx2i + dy2i)*dt,& ekh(i,j,k)*(dzh2i(k) + dx2i + dy2i)*dt ) ! if (diffnrtotl == ekh(i,j,k)*(dzh2i(k) + dxh2i(i) + dy2i)*dt) then ! imin = i ! kmin = k ! end if end do end do end do ! write(6,*) 'Peclet criterion at proc,i,k = ', myid,imin,kmin call MPI_ALLREDUCE(courtotl,courtot,1,MY_REAL,MPI_MAX,comm3d,mpierr) call MPI_ALLREDUCE(diffnrtotl,diffnrtot,1,MY_REAL,MPI_MAX,comm3d,mpierr) if (courtot <= 0) then write(6,*) 'courtot=0!' end if if (diffnrtot <= 0) then write(6,*) 'diffnrtot=0!' end if dt = min(dtmax,dt*courant/courtot,dt*diffnr/diffnrtot) timeleft=timeleft-dt dt_lim = timeleft timee = timee + dt ntimee = ntimee + 1 ntrun = ntrun + 1 else dt = dtmax ntimee = ntimee + 1 ntrun = ntrun + 1 timee = timee + dt timeleft=timeleft-dt end if end if end if end subroutine tstep_update !> Time integration is done by a third order Runge-Kutta scheme. !! !! \latexonly !! With $f^n(\phi^n)$ the right-hand side of the appropriate equation for variable !! $\phi=\{\fav{u},\fav{v},\fav{w},e^{\smfrac{1}{2}},\fav{\varphi}\}$, $\phi^{n+1}$ !! at $t+\Delta t$ is calculated in three steps: !! \begin{eqnarray} !! \phi^{*} &=&\phi^n + \frac{\Delta t}{3}f^n(\phi^n)\nonumber\\\\ !! \phi^{**} &=&\phi^{n} + \frac{\Delta t}{2}f^{*}(\phi^{*})\nonumber\\\\ !! \phi^{n+1} &=&\phi^{n} + \Delta t f^{**}(\phi^{**}), !! \end{eqnarray} !! with the asterisks denoting intermediate time steps. !! \endlatexonly !! \see Wicker and Skamarock, 2002 subroutine tstep_integrate use modglobal, only : ib,ie,jb,je,kb,ke,nsv,dt,rk3step,e12min,lmoist,timee,& iinletgen,ltempeq,idriver,BCtopm,BCtopm_pressure,BCxm_periodic,BCym_periodic, & ifixuinf,thlsrc,lchem,ierank,jerank,BCxm,BCym use modmpi, only : cmyid,myid use modfields, only : u0,um,up,v0,vm,vp,w0,wm,wp,& thl0,thlm,thlp,qt0,qtm,qtp,e120,e12m,e12p,sv0,svm,svp,& dpdxl,dgdt,thl0c use modinletdata, only: nstepreaddriver, irecydriver use modsubgriddata, only : loneeqn use modchem, only : chem use decomp_2d, only : exchange_halo_z implicit none integer i,j,k,n real rk3coef,rk3coefi rk3coef = dt / (4. - real(rk3step)) rk3coefi = 1./rk3coef if(ifixuinf==2) then dpdxl(:) = dpdxl(:) + dgdt*rk3coef ! if(ltempeq) then ! thlsrc = thlsrc + thlsrcdt*rk3coef ! end if ! write(6,*) 'dpdx = ', dpdxl(kb) end if ! Tendency zeroing and the substep-3 *m copies are fused into the integration ! nests (formerly ~14 separate whole-array passes per substep; see #330). ! The *p halos are 0 from allocation and only ever written on the boundary ! planes handled in the guarded blocks below, so no blanket zeroing is needed. ! The rk3step branch is loop-invariant (unswitched by the compiler). if (loneeqn) then do k=kb,ke do j=jb,je do i=ib,ie u0(i,j,k) = um(i,j,k) + rk3coef * up(i,j,k) v0(i,j,k) = vm(i,j,k) + rk3coef * vp(i,j,k) w0(i,j,k) = wm(i,j,k) + rk3coef * wp(i,j,k) e120(i,j,k) = e12m(i,j,k) + rk3coef * e12p(i,j,k) e120(i,j,k) = max(e12min,e120(i,j,k)) up(i,j,k) = 0. vp(i,j,k) = 0. wp(i,j,k) = 0. e12p(i,j,k) = 0. if (rk3step == 3) then um(i,j,k) = u0(i,j,k) vm(i,j,k) = v0(i,j,k) wm(i,j,k) = w0(i,j,k) e12m(i,j,k) = e120(i,j,k) else e12m(i,j,k) = max(e12min,e12m(i,j,k)) end if enddo enddo end do else do k=kb,ke do j=jb,je do i=ib,ie u0(i,j,k) = um(i,j,k) + rk3coef * up(i,j,k) v0(i,j,k) = vm(i,j,k) + rk3coef * vp(i,j,k) w0(i,j,k) = wm(i,j,k) + rk3coef * wp(i,j,k) up(i,j,k) = 0. vp(i,j,k) = 0. wp(i,j,k) = 0. if (rk3step == 3) then um(i,j,k) = u0(i,j,k) vm(i,j,k) = v0(i,j,k) wm(i,j,k) = w0(i,j,k) end if enddo enddo enddo end if ! scalars: n hoisted out of the spatial nest (n is the slowest-varying index) do n=1,nsv do k=kb,ke do j=jb,je do i=ib,ie sv0(i,j,k,n) = svm(i,j,k,n) + rk3coef * svp(i,j,k,n) svp(i,j,k,n) = 0. if (rk3step == 3) svm(i,j,k,n) = sv0(i,j,k,n) enddo enddo enddo enddo ! Chemistry is applied once per full time step, on the last RK3 substep. This ! used to be spelled `rk3coef == dt`, which holds only for rk3step==3 because ! rk3coef = dt/(4-rk3step) -- true but obscure, and an exact real comparison. if (lchem .and. rk3step == 3) then call chem end if if (ltempeq) then do k=kb,ke do j=jb,je do i=ib,ie thl0(i,j,k) = thlm(i,j,k) + rk3coef * thlp(i,j,k) thl0c(i,j,k) = thl0(i,j,k) thlp(i,j,k) = 0. if (rk3step == 3) thlm(i,j,k) = thl0(i,j,k) enddo enddo enddo end if if (lmoist) then do k=kb,ke do j=jb,je do i=ib,ie qt0(i,j,k) = qtm(i,j,k) + rk3coef * qtp(i,j,k) qtp(i,j,k) = 0. if (rk3step == 3) qtm(i,j,k) = qt0(i,j,k) enddo enddo enddo end if if ((BCxm .ne. BCxm_periodic) .and. ierank) then u0(ie+1,jb:je,kb:ke) = um(ie+1,jb:je,kb:ke) + rk3coef * up(ie+1,jb:je,kb:ke) up(ie+1,jb:je,kb:ke) = 0. if (rk3step == 3) um(ie+1,jb:je,kb:ke) = u0(ie+1,jb:je,kb:ke) end if if ((BCym .ne. BCym_periodic) .and. jerank) then v0(ib:ie,je+1,kb:ke) = vm(ib:ie,je+1,kb:ke) + rk3coef * vp(ib:ie,je+1,kb:ke) vp(ib:ie,je+1,kb:ke) = 0. if (rk3step == 3) vm(ib:ie,je+1,kb:ke) = v0(ib:ie,je+1,kb:ke) end if if (BCtopm .eq. BCtopm_pressure) then ! do i=ib,ie ! do j=jb,je ! ! w0(i,j,ke+1) = w0(i,j,ke) - dzhi(ke)*((u0(i+1,j,ke)-u0(i,j,ke))*dxfi(i) + & ! ! (v0(i,j+1,ke)-v0(i,j,ke))*dyi) ! ! if (myid ==0 .and. (i==32 .and. j==1)) write(*,*) rk3coefi*(w0(i,j,ke) - dzhi(ke)*((u0(i+1,j,ke)-u0(i,j,ke))*dxfi(i) + & ! ! (v0(i,j+1,ke)-v0(i,j,ke))*dyi) - wm(i,j,ke+1)), & ! ! dpdztop(i,j), & ! ! 2*pij(ke)*dzhi(ke+1) ! ! ! ! wp(i,j,ke+1) = rk3coefi*(w0(i,j,ke) - dzhi(ke)*((u0(i+1,j,ke)-u0(i,j,ke))*dxfi(i) + & ! ! (v0(i,j+1,ke)-v0(i,j,ke))*dyi) - wm(i,j,ke+1)) ! ! wp(i,j,ke+1) = 2*pij(ke)*dzhi(ke+1) ! end do ! end do w0(ib:ie,jb:je,ke+1) = wm(ib:ie,jb:je,ke+1) + rk3coef * wp(ib:ie,jb:je,ke+1) wp(ib:ie,jb:je,ke+1) = 0. if (rk3step == 3) wm(ib:ie,jb:je,ke+1) = w0(ib:ie,jb:je,ke+1) end if ! Write some statistics to monitoring file if ((myid==0) .and. (rk3step==3)) then open(unit=11,file='monitor'//cmyid//'.txt',position='append') if (iinletgen == 1) then write(11,3001) timee elseif (idriver == 1) then write(11, '(I4)') nstepreaddriver write(11, 3001) timee, u0(irecydriver,1,32) ! elseif (idriver == 2) then ! write(11, '(I4)') nstepreaddriver ! write(11, 3001) timee, storetdriver(nstepreaddriver), u0(irecydriver, 1, 32) else write(11,3001) timee end if 3001 format (13(6e14.6)) close(11) if (ifixuinf == 2) then open(unit=11,file='dpdx___.txt',position='append') write(11,3002) timee,dpdxl(kb) 3002 format (13(6e20.12)) close(11) if (ltempeq) then open(unit=11,file='thlsrc.txt',position='append') write(11,3002) timee,thlsrc close(11) end if end if endif end subroutine tstep_integrate end module modtstep